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SUMMARY 


Transonic Euler/Navier-Stokes computations are accomplished for wing-body flow fields using 
a computer program called Transonic Navier-Stokes (TNS). The wing-body grids are generated 
using a program called “ZONER,” which subdivides a coarse grid about a fighter-like aircraft 
configuration into smaller “zones,” which are tailored to local grid requirements. These zones can 
be either finely clustered for capture of viscous effects, or coarsely clustered for inviscid portions 
of the flow field. Different equation sets may be solved in the different zone types. This modular 
approach also affords the opportunity to modify a local region of the grid without recomputing the 
global grid. This capability speeds up the design optimization process when quick modifications 
to the geometry definition are desired. The solution algorithm embodied in TNS is implicit, 
and is capable of capturing pressure gradients associated with shocks. The algebraic turbulence 
model employed has proven adequate for viscous interactions with moderate separation. Boundary 
conditions are treated explicitly with some necessary restrictions on the solution domain. Proof of 
concept has been demonstrated with solutions for a General Dynamics F-16A-like configuration. 
Comparison of computed pressure distributions to wind tunnel test data shows good quantitative 
agreement. Comparison of computed particle traces shows qualitative agreement with both wind 
tunnel and water tunnel flow visualization of vortical phenomena. These results confirm that 
the TNS program can successfully be used to simulate transonic viscous flows about complicated 
three-dimensional geometries. 
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1. INTRODUCTION 


1.1 Objective 

For many years a goal of computational aerodynamicists has been to accurately compute 
the fluid flow around highly maneuverable aircraft configurations at transonic speeds. This goal 
required solving the necessary set of equations within an adequately fine grid of solution points 
which would correctly capture the significant physical flow phenomena. This dissertation describes 
the transonic computation of the flow around a General Dynamics F-16A derived wing-strake- 
fuselage, using the thin-layer form of the Navier-Stokes equations, and the Euler equations. The 
computed solutions accurately portray the vortex interaction off the strake with the wing flow 
region at moderate angles of attack (fig. 1), and is a significant step toward achieving realistic 
flow-field simulation for high-performance aircraft. 

With full-aircraft configuration solutions achieved, including flow through inlets and exhausts 
and external stores properly gridded, several significant improvements to the aircraft design pro- 
cess will have been accomplished. First, since access to supercomputers is expanding at a sig- 
nificant rate, with the government, industry, and the academic community greatly increasing 
investment in this hardware, more designers and potential designers can use these tools for de- 
tailed aerodynamic evaluations. The pool of adequately equipped designers has been limited to 
date. Wind tunnels have been the the sole verification recourse to analytical methods, which 
were limited to relatively simple shapes and uncomplicated flow regimes. Wind tunnels for full 
aircraft configuration evaluations, with their accompanying high expense for construction and op- 
eration, and required expertise to competently use, are well beyond the budgets and operational 
capabilities of all but the government and the largest airframe manufacturers. 

Second, the computational tool is versatile, in that it allows an easy reconfiguration of the 
aircraft being modeled, subject to the constraints and difficulties of grid generation and modifica- 
tion. Whereas major expense and time would be necessary for The reconfiguration of a comparably 
instrumented wind tunnel model, a designer may easily reconfigure a mathematical model to try 
new and exotic configurations; this can now be accomplished without fear of tunnel blockage 
problems or model structural prohibitions and without the need to lobby up the chain of com- 
mand to finance reconfiguration for a high-risk design. In an era proceeding from close-coupled 
canard applications, forward-swept, and oblique wings, a forgiving design tool allows room to 
maneuver and experiment. 

Third, because of wind tunnel wall interference, sting interference, and instrumentation lim- 
itations, wind tunnels cannot completely evaluate even some finalized concepts prior to flight 
test, such as the inboard wing flow of the X-29 forward-swept wing aircraft. Though a sting or 
wing-tip wall mount can be designed to minimize interference for a given attitude, it will always 
require a mechanism for attachment to the tunnel wall or floor, which will distort the flow field. 
Exotic programs such as magnetic suspension of models to allow sting-interference-free wind tun- 
nel testing have been discarded as impractical or exorbitantly expensive to develop, much less 
operate. Single-cable model suspension is routinely used for free-flight. tests in the NASA Langley 
Research Center 30- by 60-ft wind tunnel, but the models are expensive and the testing process 
is manpower intensive. The full aircraft computational grid solves all these interference problems 
with increasing confidence in its accuracy, and at a progressively reasonable price. 

Possibly the most significant improvement, to the design process rendered by computational 
solutions is the ease with which the data can be analysed and flow-field visualization obtained. 
The use of sophisticated software on workstations designed specifically for graphics representation 
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allows a fluids engineer with minimal graphics background the capability of easily displaying 
isolines, isosurfaces, particle traces, and a host of other information for an arbitrary solution. 


1.2 Background 

A Reynolds-averaged form of the Navier-Stokes equations, which were first introduced by 
Navier (ref. l)and Stokes (ref. 2) in the early 19th century, is employed in this study because it is 
the least complicated CFD formulation necessary to represent the flow-field physics which occur 
in the current arena of high-performance aircraft aerodynamics. The choice of what solution 
set of aerodynamic equations to solve has, and will continue to be driven by tradeoffs between 
how accurately the physical processes will be captured for a geometry of given complexity versus 
computational cost. Chapman (ref. 3) made a thorough study of these tradeoffs in which he 
categorized four major stages of approximation to the full Navier-Stokes equations: (1) linearized 
inviscid, (2) nonlinear inviscid, (3) Reynolds-averaged Navier-Stokes, and (4) full Navier-Stokes 
(with specified limitations). 

The Stage 1 approximation was the simplist and earliest employed to compute aerodynamic 
properties for complex geometries. Typically termed “paneling methods,” because the geometry 
is modeled by contiguous surface panels, this approach uses the potential wave equation. For 
comparison the full Navier-Stokes equations conserve mass, momentum and energy and use 60 
partial derivative terms in three dimensions. The potential equation, which represents a trun- 
cation of the earlier paragraph, may be derived from mass conservation with an irrotationality 
assumption, or from momentum conservation, and uses only three partial derivative terms in three 
dimensions. This approach has been successfully applied to compute subsonic flow over bodies 
without flow separation to provide pressure distributions, lift, vortex drag and side forces. It can 
also provide wave drag computation for supersonic flow over slender bodies. The utility of this ap- 
proach is demonstrated by calculation of Space Shuttle Orbiter/B747 separation studies (ref. 4), 
and F-4 store separation calculations (ref. 5). The primary limitations are the assumptions of 
irrotationality and inviscid and isentropic flow. 

The Stage 2 approximation, which in full formulation is also known as the Euler equations, 
• contains 27 of the partial derivative terms of the full Navier-Stokes equations, dropping only the 
viscous terms. This approach can successfully be used to simulate flow fields including nonisen- 
tropic processes such as strong shocks, as long as the Reynolds number is relatively high and the 
flow-field is essentially free of boundary-layer separation. For full-scale aircraft or wind tunnel 
evaluations, the Reynolds number condition is typically satisfied, but the separation limitation is 
a significant restriction on the applicability of this method. Calculations of transonic flow-fields 
around an F-16A geometry in a companion study (ref. 6) have had limited success in properly 
capturing shock phenomena. Although Euler methods will generate a vortex near a sharp edged 
surface at incidence to freestream flow (ref. 7), this occurrence is viewed as fortuitous and mislead- 
ing. This vortex disappears with grid refinement, and when present produces spurious vorticity 
and entropy (ref. 8). Since an essential aspect of F-16 high angle-of-attack maneuvering perfor- 
mance is the proper capture of strake and wing generated vortices, the Euler equations are by 
themselves inadequate. 

Chapman’s Stage 3 approximation represents the most complicated formulation of the Navier- 
Stokes equations which may be currently computed for complex geometries. The Reynolds- 
averaged Navier-Stokes equations contain all 60 partial derivative terms of the full Navier-Stokes 
equations. They are averaged over a time interval which is small compared to the mean changes in 
the flow-field, but large compared to the turbulent eddy fluctuations. Since the computational grid 
cannot, with present computers, be made fine enough to capture the smallest scales associated 
with the these turbulent eddies, models must be used to simulate the turbulent interactions 
(ref. 9). This modeling process is the Achilles heel of Stage 3 calculations. Most current solutions 
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of complicated geometries utilize a simple eddy viscosity model, typically the one developed by 
Baldwin and Lomax (ref. 10). This model has proven versatile for three-dimensional Navier- 
Stokes applications, but is limited to flows of moderate separation. Turbulence models which 
utilize a single differential equation for turbulence kinetic energy, or two differential equations 
which compute turbulence kinetic energy and another property, such as the turbulence energy 
dissipation rate, are termed one-equation and two-equation models, respectively. It is thought 
that these more expensive models compute skin friction drag more accurately, but they have not 
been widely accepted. Johnson and King (ref. 11) have developed a less expensive formulation 
based on Reynolds shear stress that may capture massive separation accurately at only a modest 
increase in computational cost above algebraic models. Successful computation of complicated 
three-dimensional flow-fields is a recent phenomena, with examples including the supercritical- 
wing flow calculations of Holst et al. (ref. 12) and Kaynak et al. (ref. 13), the elliptic missile 
computations of New’some and Adams (ref. 14), the X24-B computations of Shang and Scherr 
(ref. 15) and the airliner wing-fuselage calculations of Fuji and Obayashi (ref. 16). Shuttle- 
like flows were computed by Fuji and Kutler (ref. 17), and what were probably the first full 
aircraft configuration flow-fields computed with any form of the Navier-Stokes equations were the 
supersonic Space Shuttle calculations of Rizk and Shmuel (ref. 18). 

The Stage 4 approximation refers to solutions which directly simulate the large eddies ob- 
served to occur in the turbulent energy cascade, while smaller eddies are modeled. This approach 
is taken because the large eddies extract energy from the mean flow, and transport the momen- 
tum and energy in an anisotropic fashion unsuitable for statistical analysis. However, the small 
eddies tend to dissipate the energy they receive in an isotropic fashion, and transport little if any 
turbulent momentum or energy in a process which lends itself to modeling. Although this partial 
modeling process saves computation at the smallest turbulence scale level, the proper capture 
of the large eddies is still extremely expensive. The application of large eddy simulation (LES) 
or related techniques to Navier-Stokes calculations of realistic aircraft configurations will not oc- 
cur for many years. The survey paper of Ferziger (ref. 19) can be consulted for perspective on 
LES developments and potential. The first large eddy simulation that accurately computed the 
viscous sublayer turbulence for channel flow was by Moin et al. (ref. 20). Wray extended these 
solutions from incompressible to compressible flow (ref. 21). The possible further implementation 
of coherent struture in turbulence modeling is suggested by Coles in his Dryden Lecture oix the 
subject (ref. 22). 

One approach of capturing viscous phenomena not addressed by Chapman (ref. 3) but cer- 
tainly relevant to the design process is the use of viscous-inviscid interactive procedures. These 
procedures utilize the Boundary Layer equations, first derived by Prandtl (ref. 23) in 1904, which 
are an approximation to the Navier-Stokes equations. For high Reynolds number flows they are 
applicable in thin boundary layers, and are simpler to solve because the momentum in the di- 
rection normal to the surface can be dispensed with through order-of-magnitude analysis, and 
pressure in that direction can be treated as constant. This pressure must be specified, which is 
typically accomplished by the application of an inviscid formulation outside the boundary layer, 
with an interactive interface. This approach is much less expensive than Navier-Stokes variants, 
but it breaks down for massive separation or strong viscous-inviscid interactions. A paper by 
Steger and Van Dalsem (ref. 24) discusses the suitable applications of this approach. 


1.3 Motivation 

As stated in the initial objective, the goal of this research is to compute transonic flow 
around a General Dynamics F-16A derived wing-strake-fuselage utilizing the thin-layer form of 
the Navier-Stokes equations, and the Euler Equations. The Reynolds-averaged approach was 
taken in order to compute the separation effects which would result at the high-alpha flight 
attitudes the F-16 was designed to sustain. This separation is a controlled mechanism on the 
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strake, which is designed to cause a strong primary vortex on the inboard wing- fuselage. This 
“nonlinear” st. rake-generated lift keeps the F-16 controllable long after a conventionally designed 
wing fuselage would have stalled. Other separation which requires capture is that on the aft 
wing, which is often shock induced in transonic flight. It is also controlled in the sense that that 
supercritical wing was designed to delay the shocks which trigger that separation as far aft as 
possible. Wing tip separation is expected where the wing tip vortex causes inboard flow at high 
angle of att ack. For higher angle of attack there is also a wing generated vortex at the leading edge, 
which is once again formed by separated flow. The motivation for properly capturing these flow 
physics is both specific and general. The F-16 product line at General Dynamics has not attained 
controlled flight at the flight attitudes which its propulsion system will support. Therefore, an 
attempt is being made to reshape the forward fuselage and strake to realize controllable flight at 
higher angle of attack. The general importance of demonstrating the ability to simulate these 
physics is that all modern fighter designs, especially those being considered for the Advanced 
Tactical Fighter (ATF) competition, will be designed to use strake-, chine-, canard- or leading- 
edge-extension (LEX-) generated nonlinear lift. This solution approach is one of the few which 
can properly simulate those flows. 


1.4 Approach 

This research was carried out as part of the Transonic Navier-Stokes (TNS) project, in the 
Applied Computational Fluids Branch (RFA) at NASA Ames Research Center. The Reynolds- 
averaged Navier-Stokes equations are solved near the aircraft surface using the thin-layered as- 
sumption originally suggested by Lomax (ref. 25) and Baldwin and Lomax (ref. 10), and first 
applied by Steger (ref. 26). These equations are embodied in the ARC 3-D (ref. 27) solution ker- 
nel, which is based on the Beam- Warming implicit numerical method for compressible solutions 
(ref. 28), as modified by Pulliam and Chaussee (ref. 29) to a scalar pentadiagonalized inversion 
matrix for improved solution efficiency. The Baldwin-Lomax turbulence model (ref. 10) is used 
to calculate turbulent eddy viscosity. The Euler equations are solved in the essentially inviscid 
regions away from the aircraft.. 

The grid construction is the natural progression of a zonal approach developed by Holst et al. 
(ref. 30) for the solution of three-dimensional (3-D) wings, in which a relatively coarsely spaced 
flow field is subdivided into two types of zones: Those on the surface are treated as viscous and 
clustering is provided normal to the surface. The remaining zones away from the surface have 
coarser spacing and are treated inviscidly. This subdivision was carried out with 4 zones for 
3-D wings, but 16 or more zones are required for the wing-strake-fuselage F-16-like geometry. 
During an iteration the flow is solved in a given zone, after which boundary conditions are passed 
through adjacent faces to the zones following sequentially, as they are brought into the central 
memory and advanced one time step. These solutions were computed on a Cray XMP/48, and 
the results were graphically displayed on IRIS workstations. A CRAY 2 computer system was 
used interactively with the workstation for particle trace generation. 

The computed results are compared to wind tunnel pressure data for quantitative assessment. 
Computed surface particle traces are compared to photographed surface oil-flow visualization, and 
unrestricted particle traces are compared to photographed water tunnel dye-flow visualization for 
qualitative assessment. 
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2. GRID GENERATION 


2.1 Introduction 

There are two major reasons for employing a zonal approach to describe a complex aircraft 
geometry: The first is that this approach allows the complicated whole to be reduced to a sum 
of simple parts, each of which can be tailored for its specific requirements within the overall 
geometry. This concept includes the option to vary the nature of the equation set being solved 
in a given zone, depending on the physics expected in that region. The second advantage is the 
ability to contain the necessary variable information for all the grid points of an individual zone 
in the computer’s main memory while the solution for that zone is being accomplished. There is 
no way that a full-aircraft grid of adequate resolution to model a complex fighter geometry can 
be handled in the main memory of any but the largest of present day supercomputers. Other 
advantages possible with zonal approaches include the opportunity to modify the geometry of 
a specific zone or zones, for instance replacing one wing geometry with another, while retaining 
the original fuselage, strake, tail surfaces, etc. This would allow rapid design iteration compared 
to the process of generating a whole new grid for each component modification. Zoning also 
provides the possibility to locally refine a grid to better capture specific physics, even if the 
original component shape is retained. Furthermore, this zonal scheme has the potential to accept 
a wide variety of configurations, not just an F-16A geometry. With a relatively simple input, map 
of geometry-specific information it can generate a zoned, automatically interfaced composite grid 
ready for solution with the Transonic Navier-Stokes (TNS) Code. 

2.1.1 Zonal predecessors.- The notion of tailoring modules or zones of a solution domain 
such that specific regions of a flow field can be treated specially has honored precedent in aero- 
nautical history (ref. 31). Prandtl’s early conception of a boundary layer, and the resultant 
development of a viscous region bordered by an inviscid region within which viscous effects were 
not significant, was the the genesis of this line of thinking (ref. 23). For Prandtl this meant that 
the flow immediately adjacent to a solid surface could be treated as at rest, and all friction effects 
would be contained within a thin region near that surface; a perfect fluid assumption could be 
applied outside this region for the incompressible flow he was considering. The implication for 
computational fluid dynamicists of today is that viscous effects, which require higher levels of 
both computational cost, and solution complexity, can be constrained to thin zones near solid 
surfaces. 

As computer capabilities have expanded, the generation of smooth and accurate surface 
and flow-field grids to describe more complicated geometries has become a significant, if not 
determining, factor in the successful solution of these flow-fields. Zonal approaches to meet this 
challenge have been generated or proliferated rapidly. Boppe (ref. 32), Lee and Ruppert (ref. 33), 
and Atta and Vadyak (ref. 34) have used them to compute transonic flow-fields around realistic 
aircraft configurations using potential methods. Zonal approaches to solve the Euler equations 
have been used by Benek et al. (ref. 35), in an embedded or overset grid technique, and by 
Hessenius and Pulliam (ref. 36) and Hessenius and Rai (ref. 37) using patched grid techniques to 
solve the Euler equations for canard wing interactions. More recently, zonal approaches have been 
successfully applied to an Euler equation solution for an F-16A complete aircraft simulation by 
Karman, Steinbrenner, and Kisielewski (ref. 6). This ‘‘complex block grid system” approach has 
resulted in a smooth and well-ordered depiction of the solution volume around a highly detailed 
surface geometry, to include flow-through inlet, and exhaust, and tail surfaces. More comparisons 
of this approach and its results to the present efforts will follow. 
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2.1.2 Fundamental grid properties.- An examination of the guiding principles of grid 
generation, as demonstrated with simple two-dimensional shapes, is appropriate prior to detail- 
ing the nuances of the zonal method developed and employed herein. Solution efficiency is an 
important criteria in determining a code’s eventual usefulness. One facet of efficiency increase is 
to use a general transformation to map the flow field from the “physical,” or Cartesian domain to 
the “computational” domain (fig. 2). A general transformation is the vehicle used to straighten 
and equally space out the curved and stretched geometry of the physical domain into suitable 
uniformity for the computational domain. This is necessary because most bodies of interest in 
physical space contain curves. Physical grids also cluster solution points near areas of high curva- 
ture, both to correctly define the curved geometry, and to accurately capture the typically high 
gradients of the computed properties in those regions. The computational domain is character- 
ized by unit lengths between points and orthogonality between connecting lines. Computational 
speed is then increased because standard unweighted differencing schemes may be used in the 
solution algorithm. This uniformity of distances and intersection angles results in squares for 
two-dimensional computational space and cubes for three-dimensional computational space. 

The “computational cube,” is defined by a solution point, or node at each of its eight corners. 
There is a one-to-one mapping of point to point and cell to cell, except in regions where there are 
specially handled singularities; there one physical point may map to many computational points. 
The simplification that a zonal approach brings to this mapping process for a complicated aircraft 
geometry can be put into perspective with in analogy: Imagine stretching an entire animal hide, 
including torso and all extremities into a single square shape, compared to joining a series of 
square shaped pieces, corresponding to various body parts, attached simply at common edges. 
The sum of simple shapes is a more tractable problem both for the hide and for an aircraft surface. 

The choice of body-fitted versus non-aligned solid surface interface deserves special comment 
(fig. 3), as does the choice between patched and overlapped, or embedded zones for zonal interfaces 
(fig. 4). Body-fitted grids are used in the present zoning methodology because the termination 
of the grid as it conforms to the surface provides a simple reference area against which to cluster 
points normally. This is necessary to properly capture viscous boundary layer effects. Different 
considerations drive the choice of interface schemes at common zonal boundaries. The fluid which 
flows through these interfaces must successfully conserve physical properties, in order to correctly 
reflect transiting flow features such as shocks, slip surfaces and vortical phenomena. The two 
fundamental choices here are between patched and overlapped grids. In patched interfaces, one 
zone ends where another begins, with the only contact in a single shared line or surface. The 
alternative is to have some overlap of the grids, either complex as in the overset approach of the 
“Chimera” scheme (ref. 35), or simple, with a specified number of overlapping cells, as in the 
present approach. As depicted in figure 4 for patched and overlapped grids, the interface may be 
either continuous or discontinuous in space, function, slope or any combination of the three. 

For a two-dimensional airfoil geometry, three mesh topologies can be used to demonstrate 
the coordinate transformation from physical to computational space, and in the process highlight 
the salient features of each of these approaches. Any of the three could have been applied to 
the planes normal to the spanwise direction for the wing region of the present wing-body grid. 
The first airfoil grid, depicted in figure 5(a), demonstrates a C’-mesh, as it is wrapped around the 
airfoil both in physical space and converted to its “computational square.” Likewise figures 5(b) 
and 5(c) show an O-mesh and an H-type, or Cartesian mesh, along with their transformations for 
the same airfoil. Comparative evaluation of the approaches shows that the C’-mesh and O-mesh 
topologies allow simple control of grid refinement both normal and tangent to the surface near 
the airfoil leading edge, and avoid the singularity inherent in the H-mesh topology. At the trailing 
edge the C-mesh and the H-mesh topologies facilitate clustering normal to the surface and to a 
wake which might extend off the trailing edge, as suggested by the Kutta condition. The O-mesh, 
although useful for defining a rounded trailing edge, poorly defines the wake region. Since the 
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proper capture of viscous wake effects is one of the objectives of the TNS solution approach, an 
O-mesh would be self-defeating for the airfoil planes of the wing grid structure. 

For three-dimensional wing specification, the first letter designation will stand for the mesh 
topology in planes normal to the spanwise orientation, i.e. the above described airfoil planes. The 
second designated for planes normal to the streamwise direction. For example 0-H would infer an 
O-type grid in the airfoil plane, and an H-type grid in the planes normal to streamwise, showing 
wing cross-sections. When the fuselage mesh topology is referred to, the first letter designated 
refers to the grid orientation on the symmetry plane, while the second letter refers to the grid 
orientation in planes normal to streamwise, or fuselage cross-sectional planes. Comparison will 
show why the H-H topology was chosen for the wing grid structure, and an H-O-mesh topology 
was chosen for the fuselage grid structure. 

The difficulties of interfacing a C-0 or C-H wing-mesh topology with a fuselage grid structure 
precludes its use in simple overlapped grid schemes. The Cartesian nature of the H-H wing-mesh 
topology allows for a simple interface with the fuselage-mesh at the wing-fuselage juncture, as 
shown in figure 6 . 

Both wing-fuselage constructions have an H-mesh wing topology in this cross-sectional plane, 
but figure 6 (a) reflects an O-mesh fuselage topology, whereas figure 6 (b) depicts an H-mesh 
fuselage topology. Both approaches adequately cluster in the direction normal to the fuselage 
surface in areas of gentle curvature. However the H-mesh fuselage topology disperses points too 
coarsely spanwise in the wing-fuselage juncture. The O-mesh topology provides positive control 
of normal clustering in both the 77 (radial and spanwise) and ( (circumferential and wing normal) 
directions. This problem is not significant where wing-fuselage blending blurs the interface, as in 
the forward chordwise portion of the F-16A configuration. However, where a distinct corner does 
exist, such as near the F-16A trailing edge where the wing is attached to the fuselage “shelf,” the 
refinement problem is serious (fig. 7). 

2.1.3 Grid efficiency.- A final general area of grid discussion concerns the notion of “grid 

efficiency.” A helpful lerm for comparison purposes is the mesh efficiency ratio (MER), which 
is defined in reference 38 as “the number of grid points on the aerodynamic surface of interest 
divided by the number of grid points in an average two-dimensional surface of the same three- 
dimensional grid.” This term reflects the ability of a given mesh topology to place points on a 
surface for a given number of points in the three-dimensional flow field, and higher MERs imply 
more efficient point dispersal. The denominator of this term is simply the total number of points 
in the volume of interest taken to the two-thirds power. Inspection of table 1 ; which is taken 
from reference 38, shows that three-dimensional grids which employ O-meshes both chordwise 
and spanwise for three-dimensional wings have higher MERs than grids employing C-meshes or 
H-meshes. The C-meshes disperse many points in the wake region, clustering normal to a nominal 
wake centerline. The H-meshes are relatively wasteful, dispersing points upstream normal to the 
centerline and outboard of the wing tip in addition to wake clustering. The salvation for the H- 
mesh topology is the ability to truncate the streamwise or spanwise extent of the grid clustering 
by placing a zonal boundary at an appropriate distance from the surface of interest. It might be 
placed several cells upstream of the leading edge, outboard of a tip or downstream of the trailing 
edge for a wing, or similarly with respect to the nose and base of a fuselage. As table 1 indicates, 
this zonal correction of H-H-mesh topology makes it competitive with the other topologies in grid 
efficiency for a simple three-dimensional wing. Furthermore, if grid refinement in wake regions is 
desired, an H-H mesh allows easy extension of the necessary clustering downstream or outboard. 
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TABLE 1.- COMPARISON OF SEVERAL ISOLATED- WING MESH 
TOPOLOGIES USING THE MER EFFICIENCY PARAMETER 


Mesh topology 

MER (inviscid) 

MER (viscous) 

H-H 

0.38 

0.24 

C-H 

0.50 

0.32 

C-0 

0.79 

0.50 

O-H 

0.79 

0.50 

0-0 

1.26 

0.79 

H-H (zonal wing) 

1.13 

0.71 


The combination of H-O-mesh topology for the fuselage, and H-H-mesh topology for the wing 
resulted in excellent grid efficiency for the composite wing body. The two F-16A cases listed in 
table 2 are for the baseline 16-zone configuration and for an 18-zone configuration, in which two 
extra zones were used to cluster points circumferentially near the strake edge. The wing MERs 
and fuselage MERs displayed in the table represent MER calculations computed locally for the 
grid regions specified. As the MERs indicate, the addition of more points to the volumetric total 
lowered the 18-zone wing MER relative to the 16-zone version, but the additional points on the 
strake raised the fuselage MER, resulting in an overall improvement in grid clustering efficiency. 

TABLE 2.- COMPARISON OF F-16A WING-BODY COMPOSITE MESH 
TOPOLOGIES USING THE MER EFFICIENCY PARAMETER 


Mesh topology 


MER (viscous) surface points 


16-Zone F-16A (total points = 271,473): 
H-H (wing) 

H-0 (fuselage) 

(composite wing-body) 
18-Zone F-16A (total points = 304,134): 
H-H (wing) 

H-0 (fuselage) 

(composite wing-body) 


1.11 

4636 

1.48 

6206 

2.59 

10842 

1.01 

4636 

1.82 

8367 

2.82 

13003 


I 

2.2 Wing-Body Zonal Grid Generation 

The development of the F-16A wing-body flow-field grid, the nodes of which were used as the 
solution points by the TNS program, involved a three- step process of surface geometry definition, 
base flow-field grid development from the newly defined surface geometry, and then zoning . The 
zoning process is accomplished via a program called “ZONER," which subdivides the coarse flow- 
field grid logically into simple local regions and then clusters points approximately normal to 
no-slip surfaces. 

i 

2.2.1 Surface geometry definition.- The determination of a surface geometry definition 
for this “F-l 6A-like” configuration involved considerable engineering judgment. To begin with, 
the initial surface definition as provided by General Dynamics Corporation reflected not the 
actual airplane shown in figure 8, but rather the one-ninth scale wind tunnel model used in a 
1972 test (ref. 39) at the United States Air Force's (USAF’s) Arnold Engineering Development 
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Center (AEDC) 16- by 16-ft transonic wind tunnel (16T). The pressure data provided from that 
test was used as the primary source of experimental validation for the present computed results, 
therefore it made more sense to replicate the model surface, than that of the actual aircraft. 
This wind tunnel model was actually a YF-16A prototype version, which had the same external 
surface as the later F-16A production aircraft with the exception of a smaller horizontal tail. 
There w r ere also wind tunnel data and geometry specifications available for a one-fifteenth scale 
model of another pre-F16A prototype, as tested in 1972 at the Cornell Aeronautical Laboratory 
facilities. However, that aircraft configuration had an unreflexed strake which was atypical of 
the reflexed version finally built (fig. 9). The importance of the concave reflection area to strake 
vortex generation will be discussed in Section 4. 

The model geometry was provided in the form of (x, y, z) Cartesian coordinates, grouped on 
cross-sections of the model surface to capture the curved features of the aircraft. Mr. Thomas Ed- 
wards of NASA Ames Research Center used these points in a commercially developed computer- 
aided design/computer-aided manufacturing (CAD/CAM) software system, hosted on a VAX 
11/780 computer, to accomplish the surface fitting (ref. 40). In this process details of the model 
geometry which we did not desire to compute were modified or deleted: The most significant 
modification to the wind tunnel model geometry was that the flow-through belly inlet was faired 
over smoothly with panels, to form a no-slip surface. This was done to limit the scope of the 
initial solutions, which would involve significant complexity without the necessity of tailoring 
special zones, zonal interfaces, and boundary conditions (BC) to model the inlet. Deletions made 
at this stage included removal of the vertical tail, wing-tip missile assemblies and the smoothing 
over of surface details and gaps which were too small for the grid to resolve. 

In the second phase of surface fitting the surface points were redistributed, to group the 
computer memory-limited solution nodes where they were most needed. This was accomplished 
to more accurately define geometry, or to provide resolution in regions of the flow field where 
steep gradients in the computed physical properties were expected to occur. This surface point 
distribution forms the boundary and departure point for the flow field. Figures 9(a) and 9(b), 
taken from reference 40, show the differences between the original model geometry and the node 
distribution for the surface fitted computational grid. 

2.2.2 Base flow-field grid generation.- With the surface grid in hand, there were two pri- 
mary candidate methods for developing the base flow-field grid: One was a relatively inexpensive 
marching solution procedure provided by Edwards (ref. 41) which uses a “parabolic” differencing 
scheme to approximate derivatives in the / ] and C, directions. The other was a more expensive, iter- 
ative “elliptic” solution procedure developed by Sorenson (ref. 42), which solves a set of Poisson’s 
equations. Both procedures employed the same general strategy of developing a coarse “base” 
grid of approximately 35,000 nodes which would smoothly transition the common surface grid to 
a uniform cylinder at the far-field boundaries (figs. 11 in ref. 41 and 12 in ref. 42). This coarse 
grid would then be input to the ZONER program, which functionally subdivides it, and clusters 
points as necessary. Although the parabolic approach w r as much faster, the elliptic approach 
was selected to generate the base grid, because it provides more flexibility in controlling both 
spacing and the angles at wdiich grid lines intersect boundary surfaces. It presently takes approx- 
imately 100 sec of CRAY XMP central- processing-unit, (cpu) time to create a base grid using a 
successive-over- relaxation (SOR) procedure, but. the solution time could be reduced dramatically 
by converting to an alternating direction implicit (ADI), or other faster, scheme. Figures 7 and 
13-16, taken from reference 13, show various views of the resultant, base grid topology: figure 7 
depicts front-quarter and rear-quarter views of the surface grid. It is noteworthy that the hori- 
zontal tail surfaces were retained for this initial base grid development , but discarded prior to the 
first solution as an unnecessary complication to the initial configuration. Figures 13-15 display 
cross-sectional, planform and airfoil cuts of the base grid. Figures 16(a) and 16(b) depict the 
forebody fuselage, with special attention given to the collapse of the surface grid to a centerline 
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axis in front of the aircraft nose. The grid normal to this singularity axis is formed with concentric 
rings to fashion “nose pencils,” as depicted in figure 16(b). 

2.2.3 Zonal structure.- The program which subdivides the base grid into near-field zones 
and far-field zones is named “ZONER.” Its construction was the natural progression of a zonal 
approach developed by Holst et al. (ref. 30), Flores et al. (ref. 43), and Kaynak et al. (ref. 44) for 
the solution of three-dimensional wings. In the “wing” version of ZONER, a relatively coarse base 
grid was subdivided into four grids as shown in figure 17. The coarse, ZONE 1 extends from the 
cubical outer boundary of the flow field inwards to a finer, inviscid ZONE 2 which it surrounds, 
much like a box within a box. Because of the cubical outer-boundary shape, implementation of 
cubical wind tunnel boundary conditions in addition to free-stream outer boundary conditions was 
straightforward. There is a simple overlap of one or two grid cells between these two inviscid zones 
on common surfaces. Within the inner inviscid zone (ZONE 2) reside two viscous zones (ZONE 3 
and ZONE 4), which sit above and below the wing, and are designed to provide viscous clustering 
normal to the surface. In comparison, the wing-fuselage zoning process described below creates 
16 to 18 zones which interact with each other on single faces, with no grid totally encompassing 
another. This approach allows for simpler individual zone shapes, but requires more of them to 
define an equivalent geometry. 

2. 2 . 3.1 Specification of wing-body zonal regions: The “wing-body” version of ZONER, 
hereafter referred to simply as ZONER, follows the same philosophy of surrounding near-field 
viscous zones with far-field inviscid zones as the wing-version of the ZONER code, but does so 
without embedding any zone within another. All zones in the F-16A grid are simply connected 
volumes in the physical domain which convert to simply connected cubes in the computational 
domain. Figure 18 is a schematic which shows where boundaries lie between zones which are 
adjacent to the aircraft surface, or which lie in the plane which contains the leading and trailing 
edges of the wing. The numbering system for the zones remains fixed once established. The 
initial computations are for a symmetric configuration; the purpose is to limit, the scope of the 
initial solutions, and to halve storage requirements. Although a zoning approach lessons the main 
memory requirement, the zonal information not currently in use in the main memory still must 
be temporarily stored. We chose to compute the flow field for the right-hand side of the aircraft . 
The x-z symmetry plane is treated as a reflecting surface, although the grid extends one base cell 
beyond this plane both at the top and bottom of the fuselage. 

Some notable features of the zonal organization can be grasped by examination of figure 18; 
when viewing figure 18 and the other schematics which follow, bear in mind that there is overlap 
at the interfaces of all but. a few of the the adjacent zones, but. that an arbitrary single boundary 
is sketched to clarify the drawing. General organizational guidelines include (1) the boundaries 
of t he zonal regions are set to allow concentration of grid points in regions where they are most 
needed, keeping in mind the dimensional limitations on the size of any one zone; (2) the zones 
are placed functionally, such that thin viscous zones cover all no-slip surfaces; (3) the zones in 
the wing- fuselage juncture regions are viscous in two directions, normal to both surfaces;(4) zones 
which are viscous normal to an approximated wake centerline plane are placed behind the wing, 
and another such zone is placed outboard of the wing tip to capture the vortex which rolls up 
there’(5) zones away from the surface are treated inviscidly; and (6) the limits of each zone are 
selected to minimize the number of interfaces with adjacent zones, and the complexity of the 
interfaces that do result. For instance, where possible there will be a one to one correspondence 
of points and their locations at. the interface, resulting in direcl injection of boundary information, 
without need for interpolation. Table 3 gives the functional groupings within which all the zones 
of the 18-zone wing-body grid would fall. 
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TABLE 3.- FUNCTIONAL GROUPINGS OF ZONES CONTAINED IN 18-ZONE 
F-16A WING-BODY GRID, AS DEPICTED IN FIGURE 18’S SCHEMATIC 


Functional Grouping 

ZONE Designations 

Viscous: 


normal to fuselage 

1, 4, 5, 17, 18 

normal to wing 

8, 9 

normal to fuselage and wing 

2, 3 

normal to wake 

12, 13 

normal to fuselage and wake 

6 

Inviscid: 


flow through 

7, 10. 11, 14, 15 

slip surface 

16 


Further examination of figure 18 demonstrates the general orientation of the curvilinear co- 
ordinate system. The £ direction, which is nominally the same as the Cartesian x direction, 
increases downstream from a zero value at the streamwise location of the wing-stake juncture. 
The variable C, is oriented circumferentially around the fuselage, from bottom to top; with this 
orientation it is also approximately normal to the wing surface. The rj direction is approximately 
normal to the fuselage and parallel to the wing, increasing radially outward and spanwise. 

The cross-sectional schematics depicted in figure 19 are provided to further explain the zonal- 
interface organization. Figure 19(a) corresponds to the cut made normal to the streamwise 
direction in the aft-fuselage region, shown as “cut a” in figure 18. Figures 19(b) and 19(c) apply 
likewise to cuts through the wing-fuselage and fore-body fuselage. There is a general symmetry 
of the zonal organization around the horizontal axis, such that ZONES 6. 7, 12, 13, and 14 are 
split by the axis. ZONES 2, 4, 8, and 10 in the upper grid half are mirrored by ZONES 3, 5, 
9, and 11, respectively, in the lower half. The vertical relationship between ZONES 1,17, and 
18 are not symmetric, but could be modified to become so by the inclusion of a zone to mirror 
ZONE 18 immediately beneath it. This was not done in the present application because of the 
TNS project’s future direction. An effort is underway to replace Zone 17, which encompasses the 
fuselage underside, by several zones which model the F-16A flow-through inlet. 

Zone 16 is unique in this configuration because it is the only inviscid grid fitted to a solid 
surface. This was done because ZONE 16 simulates the exhaust cone of the aircraft, or the model 
support sting used in the wind tunnel test. The surface was treated with a slip BC, and spacing 
■was coarse in all three directions. A future grid modification will provide a flow-through exhaust 
cone, and the crudity of this modeling detail will no longer be a cause for concern. 

In figure 20 we view zonal schematics from perspectives outboard of the wing tip. These 
schematics are taken from constant i j value surfaces (designated NKC = constant), which signify 
planes in the coarse grid structure at radial increments starting from the aircraft surface. These 
surfaces define airfoils as they slice through the wing at spanwise locations. These perspectives 
give a sensation of the nonplaner nature of the grid in physical space, which must be highly 
skewed to stretch the aircraft surface shape to a uniform circular cylinder at the outer boundary. 
This viewpoint also highlights the sharp nose of the airfoil leading edge. The small leading-edge 
radius geometry proved difficult to define with the initial streamwise grid point distribution, and 
the refinement of this nose region has required considerable grid definition effort. The process of 
viewing zonal interactions in a given two-dimensional slice of the grid was made easier through the 
use of display software entitled “SLICE.” which was developed by Ms. Karen L. Gundy of Ames 
Research Center. The desire to view various facets of the flow-field domain in two-dimensional 
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slices is driven by the complex nature of the skewed grid, and its multitudinous zonal interfaces, 
which make three-dimensional grid presentation too confusing to evaluate. 

2. 2. 3. 2 Zoning control parameters: The starting and stopping points of given zones in the 

£, 77, and £ directions (fig. 18 ), which are tied to points in the base grid, are input in a directory of 
such limit points termed a “MAP.” The MAP may be easily modified, through redesignation of 
these control points, to change the extent of given zones, or to replace zones with combinations of 
other zones to add or remove components. The clustering of points normal to the no-slip surfaces 
in viscous zones is accomplished by specifying a coarse framework by means of MAP designated 
points in the base grid. Then extra points are filled in at smoothly stretched increments via an 
automatic process. 

The base grid surface station and the station two points normal to the surface form the 
framework for this process (fig. 21). Then, fineness ratio control (RAT) inputs are specified 
to establish the manner of clustering; the RAT number, which must be 2 or greater, specifies 
the number of points from those available which are positioned in the outer of the two near- 
surface cells. The balance of the normal points are smoothly clustered normal to the surface, 
such that a minimum RAT number implies the finest clustering adjacent to the surface. The 
smooth repositioning of points is automatically accomplished with a bisection iteration scheme, 
which enforces a constant stretching ratio as it operates on calculated arc lengths. The process 
of repositioning arc lengths is repeated until two conditions were met: The resultant positions of 
the viscous grid points must be both smoothly spaced throughout the two cells, and adequately 
close to the base grid locations at the three framework points were they should coincide. This 
process provides a smooth transition from a first cell dist ance from the wall of on the order of y + 
= 5 , to coarse inviscid spacing at the outer boundary for a typical finely clustered grid. 

Figures 22 and 23 demonstrate what the actual clustered grids look like, compared to their 
coarse base frameworks, as a result of this process. Figure 22 depicts the same two-dimensional 
slice normal to the streaimvise direction shown as a schematic in Figure 19 (b). The use of different 
colors for the different zones alhyws easier perception of the interfaces between zones, especially 
where the zones overlap. When there are inconsistencies in the grid such as lines crossing each 
other, which might cause negative calculated cell volumes, the use of color speeds up both the 
discovery and correction of these problems. A zoned grid can be scanned visually for grid trouble 
spots very efficiently. Figure 23 displays the same two-dimensional slice which contains the airfoil 
as the schematic of Figure 20(b). 

The bisection iteration process can be replaced with alternate stretching methods, such as a 
Newton iteration scheme, which prescribes collocation of stretched and base points at the required 
matching locations. With the arc lengths thus established, a cubic spline interpolation scheme is 
used to specify (x,v,z) coordinate locations for the resultant grid points. This clustering process 
is accomplished in one direction at a time for those zones, for instance along the wing-fuselage 
juncture, where viscous clustering was required in two directions. 

A second function of the RAT inputs is to establish a fineness ratio in those directions, within 
each grid, which do not require viscous clustering. In these instances the RAT can be specified 
as “1,” which means that the original base grid spacing is not altered, or it can be increased to 
some larger integer, such that the fineness of spacing is increased by that factor. The control 
factor which is used to designate whether the clustering will be viscous or inviscid for each zone 
is designated IVISK for the 77 direction, which is always normal to the fuselage, and IVISL for the 
c direction, which is always normal to the wing. An IVISK valued of one means viscous normal 
to the fuselage, whereas an IVISK of 0 signifies inviscid in that direction. Likewise, an 1 V 1 SL of 
1 signifies viscous normal to the wing from above; an IVISL of -1 indicates viscous normal to the 
wing from below. The IVISL of 2 is for the special case of viscous in two directions, approaching 



a wake centerline from above and below. There are no zones in this grid which required clustering 
in the streamwise, or £ direction. 

The RAT and IVIS parameters are given in table 4 for each member of an 18 -zone sample 
grid, which had relatively coarse RAT parameters specified for a laminar flow solution. The 
designations JRAT, KRAT and LRAT refer to the fineness ratio control (RAT) parameter in the 
J, K, or L index directions, respectively. Those indices correspond to the the generalized £, 77, and 
£ directions in the generalized coordinate system. Table 4 shows various combinations of viscous 
and inviscid directions for the individual zones. 


TABLE 4 .- CONTROL PARAMETERS FOR AN 18 -ZONE GRID 


Grid no. 

JRAT 

KRAT 

LRAT 

IYTSK 

IVISL 

1 

4 

6 

1 

1 

0 

2 

4 

6 

6 

1 

1 

3 

4 

6 

6 

1 

-1 

4 

1 

6 

2 

1 

0 

5 

1 

6 

2 

1 

0 

6 

2 

6 

3 

1 

2 

7 

2 

2 

1 

0 


8 

4 

2 

6 

0 

1 

9 

4 

2 

6 

0 

-1 

10 

1 

2 

1 

0 


11 

1 

2 

1 

0 


12 

2 

2 

3 

0 

2 

13 

2 

2 

3 

0 

2 

14 

1 

1 

1 

0 


15 

2 

1 

1 

0 


16 

1 

1 

1 

0 


17 

4 

6 

1 

1 


18 

4 

6 

6 

1 

1 


2. 2. 3. 3 Dimensions and global grid characteristics: The dimensions and total number 
of points for each grid of the 18 -zone wing body are given in table 5 . There was a limitation of 
30,000 grid points available for a given grid, based on program structure. 
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TABLE 5.- ZONAL DIMENSIONS AND TOTAL POINTS PER ZONE 
FOR AN 18-ZONE GRID CONTAINING 304,134 POINTS 


Grid no. 

NJ 

NK 

NL 

Product 

1 

69 

19 

12 

15732 

2 

77 

19 

19 

27797 

3 

77 

19 

19 

27797 

4 

37 

19 

23 

16169 

5 

37 

19 

23 

16169 

6 

37 

19 

19 

13357 

<"T 

i 

35 

29 

25 

25375 

8 

77 

19 

19 

27797 

9 

77 

19 

19 

27797 

10 

37 

29 

12 

12876 

11 

37 

29 

12 

12876 

12 

37 

19 

19 

13357 

13 

73 

13 

19 

18031 

14 

53 

6 

25 

7950 

15 

9 

20 

25 

4500 

16 

6 

20 

25 

3000 

17 

69 

19 

13 

17043 

18 

69 

19 

19 

24909 


For comparison purposes the total number of points employed in the four-zone-wing-alone 
grids of references 12-14 was 149,000, and approximately 270,000 points were used for the 16- 
zone version of the F-16A wing-body. 


2.3 Zonal-Interface Scheme 

The zonal-interface scheme developed for the wing-body version of ZONER represents a signif- 
icant. improvement over the version which was used in the isolated-wing version. In both versions 
the objective was to accurately pass information from neighboring zones to the boundaries of the 
zone being solved for, with sufficient accuracy to guarantee conservation. If mass, momentum 
and energy are not. conserved through the transition, then phenomena such as shocks, which cross 
the interface, will not. be properly captured. The improvement which was made concerned the 
way the solution was interpolated for passage through discontinuous interfaces. 

The wing- body ZONER retains the simple overlap approach of its predecessor, such that there 
is a one-base grid-cell overlap between the boundaries of two adjacent zones, irrespective of what 
the R ATs may be for each zone in the direction normal to the interface. Figure 24(a), taken from 
reference 44, depicts this process between a relatively fine grid and a coarse grid to illustrate the 
point. Fine “ZONE 2,” say of RAT 2 in the £ direction, overlaps coarse “ZONE 1” of RAT 1 by 
one base grid cell, and ZONE 1 overlaps ZONE 2 in the same fashion. In this instance the fine 
zone also has higher fineness ratios in the ( and rj directions, but that need not have been the case. 
Figure 24(b) shows the cross section of the interface plane, and it is obvious that since ZONE 
2 has a “+” point at. every point that ZONE 1 has a “0” point, data transfer from ZONE 2 to 
ZONE 1 can lie achieved by simple injection at those common points. The transfer in the other 
direction is more complicated though; solutions must be interpolated from between the coarsely 
spaced grid points of ZONE 1 to provide the needed information at all the finely spaced points 
of the ZONE 2 interface plane. 
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In the wing ZONER version this was accomplished via a simple linear interpolation process 
in one coordinate direction, which in that instance was normal to the wing. In this wing-body 
ZONER, we use cubic spline interpolation, and use arc length based on the two grid directions 
normal to the interface direction as the independent variable. The process is accomplished lo- 
cally for every nodal location in the interface plane, and the approach is more accurate than its 
predecessor for two primary reasons. The higher order of spline interpolation properly captures 
the nonlinear curvature of the grid, which is highly skewed in some regions. Secondly, the use 
of arc length allows the process to be applicable for interfaces of irregular surfaces, in which 
the points in any given direction might not behave in a monotonic fashion, as required by the 
single direction approach. Since arc length is by definition a monotonic concept, the grid can 
turn in any direction, and contain inflection points, without causing interpolation inadequacies. 
This approach is thus capable of interpolating for grids of arbitrary shape with a high degree of 
confidence. 


2.4 Data Management 

Once the grid is generated and properly subdivided by ZONER, it can be accessed by the TNS 
code to start the solution. The iteration procedure begins with Zone 1, in the fuselage forebody, 
and continues sequentially through the remaining blocks, before it starts the next iteration again 
in Zone 1. Thus only one iteration is accomplished in each block per configuration iteration, and 
the rotational order through the zones does not vary. While a given zone is being solved, the 
solution information for that zone only resides in main memory. Similar information for the other 
zones is maintained in extended storage on a CRAY Solid Storage Device (SSD). Thus the main 
memory need hold only the solution and transformation Jacobian array (Q), metric quantities 
and turbulence model property arrays for one zone at a time. The SSD is functionally used in 
the same fashion as standard rotating disk extended storage. The SSD is physically made up of 
semiconductor memory, however, and as such can be accessed much more rapidly. To allow more 
space for solution points in main memory, we retain the metric quantities for only two-dimensional 
planes as needed; since these metric planes are continuously being shuffled in and out, if we had 
to rely on standard rotating disk or tape storage, we would be Input/Output (l/0)-bound, and 
the solution process would be unpractically slow. 

The SSD storage capacity has limitations though, which must, be considered. The XMP/22 
was the machine on which the TNS code was originally designed to run, so the storage process 
was allocated with its storage potential in mind. The SSD associated with the CRAY XMP/22 
had 16 million 64-bit words of memory available; that capacity could be doubled to 32 million 
words if 32-bit precision was used. This turned out to be perfectly acceptable when used for the 
metric arrays, as determined in a test case. The significance of this substitution was caused by 
the use of an Alternating Direction Implicit (ADI) algorithm as the basis for the solution method. 
Sweeps were performed in all three directions, requiring the metric quantities to be accessed in 
three different orientations. This exorbitant storage requirement was easily met by the XMP/22 
SSD in the 32-bit precision mode. 

When the CRAY XMP/22 was upgraded to an XMP/48, the same approach of storing the 
three-dimensional metric arrays on SSD, and accessing two-dimensional slices of the metrics was 
retained, even though the whole arrays could have been brought into the main memory which 
now comprised four million 64-bit words. This was done to prevent interference with other jobs, 
since XMPs are time-sharing systems. 

With the access available to the CRAY 2 system at Ames Research Center, we have the 
potential to run with all necessary information available in a main memory which can hold 256 
million 64-bit words; the CRAY 2 does not have an SSD. Since the SSD I/O transfer caused 
approximately a 30% efficiency penalty when compared to main memory access, this will cause a 
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significant drop in the I/O requirements. More significant will be the ability to refine the surface 
geometry to define smaller, though aerodynamically significant, aspects of the “real” aircraft. It 
has been suggested that this larger memory capacity could make zonal approaches superfluous, 
in terms of their utility for minimizing memory requirements. One need only consider that 
present-day full-aircraft grids can achieve resolution no finer, on average, than a credit-card-sized 
rectangle, to realize that orders of magnitude memory increase will be needed to properly model 
aircraft gun ports, sharp strake edges, and external store details. The desire to incorporate 
deterministic turbulence computations, rather than modeled turbulence data, could easily absorb 
all storage increases on the horizon. 



3. THE PHYSICAL EQUATIONS AND 
THE NUMERICAL METHOD 


3.1 Introduction 

The strong conservation-law forms of the combined Euler and Reynolds-averaged Navier- 
Stokes Equations are used to solve for the flow fields in this study. The Reynolds-averaged 
equations are simplified by using the standard thin-layer approximation for the viscous terms, as 
first suggested by Baldwin and Lomax (ref. 10), and applied by Pulliam and Steger (ref. 26). Once 
the zonal approach has been established, a wide variety of standard Navier-Stokes integration 
algorithms can be used. In this study, the basic governing equations and numerical algorithm 
including the turbulence model have been taken from Pulliams ARC3D computer code (ref. 27), 
which was developed from the Beam and Warming implicit approximate factorization scheme 
for finite-difference solutions (ref. 28). That scheme was made more efficient and robust by 
matrix modifications introduced by Pulliam and Chaussee (ref. 29), to transform block tridiagonal 
solution matrices to pentadiagonal scalars, which incorporate an improved dissipation operator. 
Boundary condition specifics and the artificial dissipation and turbulence models employed will 
also be discussed. 


3.2 The Governing Equations 

3.2.1 The Reynolds-averaged Navier-Stokes equations.- The Navier-Stokes equa- 
tions have the capacity to describe the flow of an unsteady, viscous, compressible, continuous 
fluid in the absence of body forces and electromagnetic fields. The full equation set cannot yet be 
solved for high Reynolds number flows around complicated geometries. It is too computationally 
expensive to make these grids spatially fine enough to capture the small scales associated with 
turbulent eddies; therefore, turbulence models must be employed. The implementation of these 
models typically begins with the time averaging of the Navier-Stokes equations, a process in which 
the time-dependant, quantities a're divided into a slowly varying or mean-flow part, and a high- 
frequency part. The resulting equations are averaged over a time period which is large compared 
to the fluctuating period of the turbulent eddies, but small compared to the time elapsed during 
mean changes in the flow field, as first suggested by Osborne Reynolds in 1883 (ref. 45). 

Conservation implies that mass, momentum and energy are conserved in all cells, including 
those which transit regions of extreme gradients such as shocks; this control volume concept 
extends to the total grid, which is the sum of all such small cells. These equations are solved in 
strong conservation-law form (ref. 46), which means that neither the functions of the flow variables 
nor their metrics occur as coefficients in the differential equations. This form is used in most 
instances when shock capturing is an important consideration. If the metrics were to appear as 
coefficients, the less desirable weak conservation-law form is described. The strong conservation- 
law equations thus specified can be written in Cartesian coordinates and nondimensional form 
(ref. 47): 

OQ dE dF 00 _ 3E r 
+ + ~dy + 0z~ ~dx 
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(3.3) 


The Cartesian velocity components u, v, and m are nondimensionalized by a oo (the free- 
stream speed of sound), density p is nondimensionalized by /><*,, and the total energy per unit 
volume e is nondimensionalized by p^a'^. Pressure can be found from the ideal gas law as: 


p = (7 — l)[e — 0.5p (» 2 + r 2 + ir 2 )] (3.4) 

The ratio of the specific heats, 7 , is set equal to 1.4. Also, k is the coefficient of thermal conduc- 
tivity, p is the dynamic viscosity, and A , a coefficient of viscosity from the Stokes’ hypothesis, is 
— 2/3jU . The Reynolds number is Re and the Prandtl number is Pr. 

The transformation of these equations to the computational domain (f,? 7 , () is accomplished 
in a manner which preserves their strong conservation form. The time-dependent metrics which 
occur in the most general form of the equations have been omitted for the present steady-state 
solutions, although they are presented in equation 3.7 for completeness. If the grid is allowed 
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to vary with time in an adaptive manner, these time-dependent metrics would be used. The 
transformed curvilinear coordinates are defined as: 


r = t 

t = Z{x,y,z,t) 
r] = r](x,y,z,i) 

C = 

For steady-state solutions At can be treated as a spatially varying quantity to 
solution convergence rates. The transformed governing equations are given by: 
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U, V', and W are contravariant velocity components. The viscous flux terms are 
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where the components of the shear-stress tensor and heat-flux vector in nondimensional form 
were given in equation (3.3). Here, the Cartesian derivatives are expanded in £, y, ( space via 
chain-rule relations such as: 

U x = ZxU( + r] x u v + Cr«c 

The metric terms are obtained from chain-rule expansion of x y v , etc., and solved for £*, 
etc., to give: 
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and the transformation Jacobian J is defined by 


J 1 = x^y v z < + x ( i ^z v + XT,y<:t ~ *m z v ~ - x cVv z i 

Different, approaches are needed for the Jacobian calculation in the wing-fuselage corner 
regions, because of the excessive aspect ratio of the cells which are viscously clustered in two 
directions. The Jacobians calculated from the finite-differenced derivatives in equation 3.10 are 
inaccurate where the long thin cells are highly skewed. A scheme by Chaderjian (ref. 48) which 
determines the Jacobian from consistent metric differences is both accurate and efficient for 
interior points, but cannot be used for boundary cells. Here a volumetric approach, based on 
the realization the the Jacobian is the inverse of the cell volume, is applied. This approach 
uses a centroid pyramid formula in an accurate, robust, and inexpensive manner. An alternate 
volumetric approach using trilinear quadrature was discarded as sufficiently accurate, but too 
computationally expensive. 
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3.2.2 Thin-layer approximation.- The thin-layer approximation, as first used by Pulliam 

and Steger (ref. 26), takes advantage of the observation that, for high Reynolds number flows the 
viscous effects are confined to a thin layer near rigid boundaries. The simplifying feature of this 
observation is that flow gradients vary rapidly only in the direction normal to the surface within 
this thin layer. Since extremely fine grid spacing is required to capture high gradients, it is doubly 
fortunate that it need be employed only in a thin region, and only in one direction. This allows 
us to simplify the Reynolds-averaged Navier-Stokes equations by assuming that viscous gradient 
terms in the two general curvilinear directions tangent to the solid surfaces are very small; they 
may be dropped from the formulation. 

For the wing-body grids described in this study, // is the direction normal to the fuselage, 
and ( is the direction normal to the wing surface. Clustering is required in both directions for the 
wing-fuselage junctures, as described in Chapter 2. Considering application of this assumption 
for flow in the junctures equation (3.6) simplifies to: 
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This approximation is valid only for high Reynolds number flows where the boundary layer is thin. 
Very large turbulent eddy viscosities, associated with massive flow separation, may invalidate the 
model. 


3.2.3 The Euler equations.- The Euler equations are an inviscid subset of the Navier- 
Stokes equations, and as such are also a subset of the thin-layer version, which alters only viscous 
components. They are recovered from the thin-layer formulation by dropping all viscous terms, 
which in this instance means setting the right hand-side of equation (3.11) equal to zero. 
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3.3 Numerical Method 


A fundamental algorithm decision is the choice between explicit and implicit methods: Ex- 
plicit algorithms advance from one time step to the next using information available from the 
previous time step, or earlier ones. Implicit methods, on the other hand, include information 
from adjacent points in the flow field at the same time step as the point being solved for. This 
necessitates the simultaneous solution of sets of linear equations, and results in expensive matrix 
inversions. The simplicity and cost advantages of the explicit approach are offset by limitations 
in the allowable size of the time step. If these are exceeded instability, which is a reflection of 
the system’s inability to return after a perturbation to a state of equilibrium, renders the method 
ineffective. Large time steps are desired to reach a steady-state solution rapidly, but classical 
linear stability analysis (ref. 28) predicts only conditional stability for explicit methods, such that 
small time steps must be taken. A related measure of merit is stiffness, which signifies the relative 
restriction to the time step limit imposed by a stability bound versus an accuracy bound. It is 
the wide range of length scales from fine to coarse grid regions which introduces the stiffness 
that causes stability limits for explicit methods. Similar analysis predicts unconditional stability 
for implicit methods, which makes them attractive for problems which require fine grid spacing 
for numerical resolution. Since the equations being solved are nonlinear, there are time step 
limitations even for implicit methods, but they are not as restrictive as for explicit approaches. 

3.3.1 Beam- Warming block ADI algorithm.- The finite-difference method employed in 
the ARC3D code, which is the core from which the TNS program is fashioned, is the implicit ap- 
proximate factorization algorithm in delta form, by Beam and Warming (ref. 28). This alternating 
direction implict (ADI) algorithm sweeps in each of the three transformed coordinate directions 
sequentially. Its significant improvement over previous ADI methods lay in its successful retention 
of conservation properties through proper linearization of terms. Instead of solving a sparse, but 
expensive block matrix of order of the £ dimension times the rj dimension times the ( dimension, 
three one-dimensional matrix operators can be factored out and inverted individually, resulting 
in a significant comput ational savings. The cross terms of the original operator are second-order 
accurate in time, so they can be disregarded if the chosen scheme is second-order accurate or 
less. The resultant operators, which are block tridiagonal because of the second-order central 
differencing employed, can be sequentially solved using lower-upper decomposition (LUD). 

When the finite-difference algorithm is applied to the thin-layer approximation, equa- 
tion (3.11) reduces to: 


[I + h8^A n + D< 2) ]x 

[/ + hS„B n - hRe~ 1 6 v J~ 1 N n J + D* 2) ]x 
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where 6 is the central- difference operator and A and V are forward and backward-difference 
operators; e.g., 
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/ *} \ / 2 \ 

The time step h = At corresponds to first-order time-accurate Euler Implicit and D y , D v 
and Z)^ 2) are the implicit and D ! 4) is the explicit smoothing operator. The Jacobian matrices 

A n , B n and C n are obtained in the time linearizations of E n , F n and G n . These flux jacobians 
and the viscous coefficient matrix M, which comes from the time linearization of the viscous 

vector S n+1 , are documented in appendix A. The other variables in equation 3.13 are described 
in reference 44. 


This algorithm is first- or second-order accurate in time and second- or fourth-order accurate 
in space. The stability and accuracy of the numerical algorithm is described by Beam and 
Warming (ref. 28). According to the linear analysis, the numerical scheme is unconditionally 
stable in two dimensions, but in actual practice time step limits are encountered because of the 
nonlinear nature of the equations. In three dimensions the algorithm is unstable unless artificial 
dissipation is added. 

3.3.2 The diagonal ADI modification.- Approximate factorizat.on techniques alone are 
not sufficient to make block implicit algorithms efficiently solvable for grids as large and com- 
plicated as the present requirements dictate. The obstacle of the expensive inversion of block 
tridiagonal matrices has been circumvented by Pulliam and Chaussee (ref. 29), who took advan- 
tage of the eigenstructure to diagonalize the Jacobian matrices, as suggested by Warming, Beam, 
and Hyett (ref. 49). This process converted the block tridiagonal matrix to a scalar tridiagonal 
variety. 

The Jacobian matrices, .4, B, and C are reduced through similarity transformations to: 


A = T(X ( Tr', 


B=T,k,T-\ C = W, 


-1 


(3.15) 


To accomplish this the eigenvector, or similarity transformation matrices TV, T v , and 
are factored outside the spatial derivative terms, and new combinations of implicit smoothing 
operators (Di |^, D{ \ v and Di\^) which easily diagonalize in scalar fashion are incorporated. The 

variables A*, A,,, and A^ are diagonal eigenvalue matrices. As a result equation (3.13) is trans- 
formed to tlie scalar pentadiagonal form, which reads: 


T(_(I + h6{A.£ — hDi\^)N (I + hS^ — /?ZJ,;|^)x 

P{I + hS c A c - hD i \ < )£\Q n = R n 


where relations exist between 7^, T v , and of the form 

iV = T- 1 T 7/ , , A -1 = P = T; 1 T 0 ,P~ l =Ti' i T 1) (3.17) 

The similarity transformation matrices 7^, T v , T^, and their inverse matrices are given in ap- 
pendix B. 

Since the modifications are restricted to the implicit, or left-hand side of equation (3.16), 
the converged steady-state solution will be identical to the one obtained from the unmodified 
algorithm. Further, two-dimensional, linear stability analysis shows that the same unconditional 
stability is retained. This modification reduces the scheme to first-order accurate in time, and 
adds a nonconservative feature which degrades shock capture for the time-accurate mode, but for 
steady-state solutions accuracy is not impaired. 


Flores applied scalar diagonalization to the block tridiagonal TNS version of ARC3D (ref. 50). 
He realized as much as a 40 fold decrease in CPU time requirement from the inefficiently coded 
solutions which the “wing” TNS code version had been providing (fig. 25). During this imple- 
mentation Flores discovered that specifying different Courant-Friedrichs-Levy (CFL) numbers in 
different zones significantly enhanced convergence, a technique which has since been employed 
routinely for wing-body solutions. 

Pulliam and Steger (ref. 26) recommend space varying At to improve convergence for steady- 
state solutions, especially for grids with spacings which vary from very fine to very coarse, as in 
the present case. Two spatially varying At schemes were tested to improve convergence rate in 
the wing-alone study: 

n ft 

At = -^-=- (3.18a) 

^max 

and 

At = — % (3.186) 

1 + sfj 

where J is the Jacobian of the coordinate transformation, \ max is the maximum eigenvalue of 

the flux jacobians A, B, and 6', and A i re j is a reference time step. The first time stepping 
technique performed marginally, but the latter scheme has significantly increased convergence 
rates (ref. 44). 

An alternate approach to correcting the efficiency problems of the block tridiagonal inver- 
sions is a matrix reduction method by Barth and Steger (ref. 51). This technique reduces the 
original large block tridiagonal matrix to two smaller scalar tridiagonals and one block two-by- 
two tridiagonal matrix. A significant reduction in the number of arithmetic operations needed to 
solve the approximate factorization solution is again demonstrated. 

3.4 Boundary Conditions 

For complicated wing-body grids the selection, implementation and assessment of BCs are 
as significant and challenging as any other aspect of the problem. Among the requirements that 
the chosen BCs must satisfy are that the physics of the flow be properly defined, that they be 
posed in proper mathematical perspective for the problem, and that the numerical approximation 
must be accurate and efficient without degrading the stability of the interior solution. In addition 
BCs must be averaged or otherwise cleverly specified in grid regions where singularities exist. 
Sometimes numerical BCs must be fabricated to satisfy computational stencils which extend 
beyond physical boundaries. 

3.4.1 General Specifications.- All BCs are implemented explicitly, which would seem at 
variance with the implicit flow solver. This is done because incorporation of implicit BCs requires 
modification of the inversion matrices, which is both complicated and expensive. By treating 
them explicitly, boundary conditions can be modified or changed in a modular fashion, without 
affecting the basic properties of the code solver. Previous usage has shown that this explicit 
treatment, which gives zeroth-order time extrapolation, and either zeroth- or first-order space 
extrapolation, has caused no serious discrepancies in the solution (ref. 27). 

General specifications to this wing-body geometry, which were briefly mentioned in chapter 
2, include the enforcement of no-slip conditions everywhere on the aircraft surface by setting the 
surface velocities to zero. The sting which extends aft of the aircraft tail (fig. 14) is treated 
invisicdly, and flow tangency conditions are imposed. Pressures on the surface are determined by 
applying the boundary-layer approximation, where dp/dn =; 0 , such that Pi = p> Pi is the static 
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pressure on the body surface, and p 2 is the static pressure on the first computational surface off 
the body. The surface density is found from an adiabatic wall assumption which yields pi = p 2 . 
Total energy is computed using equation (3.4). Free-stream boundary conditions are set at all 
farfield boundaries. This approach can cause perturbations in the solution if the boundaries are 
close to regions of high gradient. When this is a factor, a more robust characteristic approach 
(ref. 27) can be employed. 

3.4.2 Grid Anomaly Regions.- In the present approach, the conditions at the outflow 
boundary are found by a zeroth-order extrapolation from the last plane such that Qjmax = 

Qjmax- j • The symmetry plane BCs are more complicated: Here a zeroth-order space extrapola- 
tion is used for the density. A first-order extrapolation is employed for the x-direction Cartesian 
velocity component u and z-direction Cartesian velocity component w , while the spanwise y- 
direction Cartesian velocity component v is set equal to zero to force symmetry. Pressure is also 
found by a first-order extrapolation and the total energy is computed using equation (3.4). 

Along the extension of the fuselage forebody grid forward into the flow field, a singularity 
region is created because the circumferentially distributed points collapse onto each other, to 
form what might be viewed as a fuselage centerline axis (fig. 16). The solution values here were 
determined by taking the average value of the superimposed points. This averaged value was 
then specified for these points, as required for future computation. 

The “interface” of zonal boundaries, as described in section 2, is applied explicitly. This 
treatment causes no problems unless one attempts to run the code at lower subsonic speeds; if 
this is attempted the solutions on either side of the interface will not converge. 

There is a special treatment interface region in the £ — rj plane which splits the top half of the 
wing flow field from the bottom half (fig. 26). The zones facing each other across this plane do 
not overlap, and averaging is employed to establish the values on the plane. Where the spacing 
to the first node on either side of the plane is equidistant, this process is straightforward. Where 
the distances are not equal, as with the interface situation between viscously clustered Zone 18 
and inviscidly clustered Zone 17, a linear weighting function is employed to increase the influence 
of the solution from the nearer node. The cross section of the forebody in figure 27, and the 
enlarged detail in figure 28 display this disparity in node distances. 


3.5 Artificial Dissipation Treatment 

The solution of the Euler and Navier-Stokes equations by numerical techniques requires the 
addition of some form of artificially introduced dissipation. This is especially true for the Euler 
equations because the physical processes being modeled do not provide any natural dissipation. 
For practical purposes it is equally necessary for the Navier-Stokes equations, to dampen high 
frequencies which are caused by shocks and other nonlinearities in the flow fields, and to control 
odd/even decoupling of grid points which is typical of central differencing. The principal challenge 
to successfully employing atrificial dissipation is to apply enough in the appropriate areas to 
prevent unacceptable instabilities in the numerics, without washing out the crisp capture of the 
physical processes. 

There are two basic ways to introduce artificial dissipation: One is through the use of upwind 
differencing schemes, such as the “monotone”, total-variation-diminishing (TVD), “flux-split.” and 
“flux difference ” methods, which automatically include a dissipation mechanism. The work of 
Steger and Warming (ref. 52), Roe (ref. 53), Van Leer (ref. 54), Osher and Chakravarthy (ref. 55) 
are of this genus. These approaches were generated under the assumptions of characteristic t heory 
and wave propagation, and they are all fundamentally equivalent to central differencing schemes 
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to which dissipation has been added. With these approaches the dissipation occurs automatically, 
so the user has limited control of both the amount applied and its effect on the solution. 

The alternate approach, to use central differencing and consciously add artificial dissipation, 
has been widely used in MacCormack’s explicit scheme (ref. 56), the explicit Runge-Kutta meth- 
ods of Jameson, Schmidt and Turkel (ref. 57), the work of Rizzi and Eriksson (ref. 58), and in 
this study. When the block tridiagonal solution approach is used, an explicit fourth-difference 
dissipation model must, be employed because a fourth-difference implicit variety would alter the 
matrix structure inappropriately. A second-difference dissipation operator is then added to relax 
the stability bound introduced by the explicit operator. The use of the Pulliam and Chaussee 
(ref. 29) variation eliminates this difficulty, because the implicit fourth-difference dissipation op- 
erator diagonals may be conviently added to the converted scalar tridiagonal solution matrix to 
form the efficiently configured scalar pentadiagonal result. The stability and convergence of the 
scheme are now enhanced at a minimal increase in cost. Implicit fourth-difference dissipation can 
still be improved upon, since the solution is a nonlinear process, in regions of high gradients. Near 
shocks an implicit second-difference operator should be added to take advantage of its narrow 
three point stencil. A toggle which activates this operator locally is typically used. The combined 
difference operator is given by: 


Oils = Vf(<r J+ 1 .*,,J- + Y W + 


J 4) A 


<v { A { p 


(3.19) 


where ( and ( are the coefficients for the second- and fourth- difference operators respec- 
tively. The terms for D,[ n and D,\( have an equivalent format. 

We have used only the fourth-difference part of the operator to date. The streamwise cluster- 
ing has not yet been fine enough in the vicinity of the shocks for the five point fourth-difference 
stencil, which spans it, to cause overshoot or undershoot stability problems. Pulliam’s excel- 
lent paper (ref. 59) on the subject contains more detail on the fine points of applying artificial 
dissipation. 


3.6 Turbulence Models 

The Baldwin- Lomax turbulence model (ref. 10) has been used throughout this study, as well 
as predecessor studies, to determine a turbulent “eddy viscosity” with which to combine the 
laminar eddy viscosity in the Reynolds-averaged Navier-Stokes computations. This phenomeno- 
logical approach to the measure of turbulence, and its effect on high Reynolds number solutions, 
is problematical, in that eddy viscosity is essentially a contrived quantity which models physical 
processes it cannot hope to concisely represent. 

3.6.1 Development of turbulence measures.- The phenomenological modeling concept 
results from Boussinesque’s replacement of the molecular viscosity coefficient in the Navier-Stokes 
equations by a turbulent viscosity coefficient (eddy viscosity) in an 1877 presentation (ref. 60). 
This new quantity represented a blend of flow properties and material properties, and Prandtl’s 
(ref. 61) subsequent creation of the “mixing length” in 1925 popularized a way of viewing turbu- 
lence which has been widely used. It has, however, suppressed the development of more fitting 
alternatives. The structural approach to turbulent property calculation, as described by Chap- 
man and Tobak (ref. 62), seeks to disprove the statistical assumptions of the phenomenological 
approach that all turbulent motions are random, and therefore must be averaged in a nonde- 
terministic fashion. The adherents to this approach have used experimental observation and 
numerical calculation to quantify turbulent structures, such as the turbulent spot first reported 
by Emmons (ref. 63). But although large eddy simulation and other approaches have achieved 
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some quantification, as yet no theory has resulted which could replace the phenomenological 
approach in a broad enough way to be meaningful for practical applications. 

3.6.2 Phenomenological model candidates.- The algebraic models, also called zero- 

equation models since they require the solution of no partial differential equations, include the 
Baldwin-Lomax model and the Cebici-Smith model (ref. 64) from which the former was derived. 
The major difference between the two is that the Cebici-Smith model requires the specification of 
the BL thickness, whereas the Baldwin-Lomax model determines it from the calculated vorticity. 
The Cebici-Smith model performs well for solutions to the BL equations, with which it is eas- 
ily applied. The Baldwin-Lomax model is more generally applicable, and has been the preferred 
method for Navier-Stokes solutions. Both methods are of questionable accuracy for the prediction 
of moderately to massively separated flow regions. 

A half-equation model has been proposed by Johnson and King (ref. 11), which is so des- 
ignated because this approach uses only a single ordinary differential equation, which measures 
Reynolds shear stress. This formulation has exceeded algebraic model results for two-dimensional 
flow calculations with moderate separation, and plans have been set to modify it for use in a 
separated three-dimensional calculation. 

The two equation models, such as the k-c model, require the solution of two partial differential 
equations, which calculate the kinetic and internal energy. They are relatively expensive to apply, 
especially for fine grids. These approaches have not consistently improved the capture of turbulent 
flow features commensurate with their increased expense and, so, have found limited popularity 
in computational applications. 

3.6.3 Baldwin-Lomax model specifications.- The Baldwin-Lomax model is a two-layer 
algebraic model in which p? is given by 


pr = 


inner ? 
(^T )outcr i 


y ^ l lcrosscn<er 
y ^ Vcrossoi'er 


(3.20) 


where y is the normal distance from the wall and y C rossover is the smallest value of y at which 
values from the inner and outer formulas are equal. The eddy viscosity coefficient in the inner 
layer is based on the Prandtl mixing-length theory 


(/*T)inner = I (3.21) 

The parameter £ is the mixing length corrected with the Van Driest damping factor to account 
for the laminar sublayer 

£ = ky[ 1 - exp(-y + /A + )j (3.22) 

wdiere k = 0.4, A" 1 " — 26, and | u> | is the magnitude of the vorticity given by: 


|<x> I = \J{u y - v x ) 2 + ( C; - Wy ) 2 + ( W x - U : ) 2 (3.23) 

This vorticity calculation is typically made with difference approximations which assume an or- 
thogonal grid. If the grid is not orthogonal near a surface. Barth (in private communication) has 
suggested using an alternate approach which makes use of contravariant velocities, to improve 
the accuracy of the calculation. 
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The y + values are obtained from 


+ PwU T y 

y = 


\J PwT w y 
Vw 


(3.24) 


The eddy viscosity coefficient in the outer layer is based on the distribution of vorticity which is 
used to determine the length scale and is given by: 

(PT)outer = K Cc p P Fw ak e Fk le b(v) (3.25) 


w’here I\ = 0.0168 is the Clauser constant and Ccp = 1.6 is an additional constant. The function 
Fwake is defined by: 

VmaxF max 


Fwake = min 


or 


@WK Vmax^\)ip / F n 


(3.26) 


where Cwk — 0.25 and 


udif = (y/u* + v 2 + w 2 ) max - ( %/ u 2 + v 2 + w 2 ) min 


(3.27) 


In equation (3.27) the second term is taken to be zero (except in wakes). The quantities 
Vmax and Fmax are determined from the function 


F(y) = y \*> |[1 - exp(-y+ /A + )] (3.28) 

The quantity F max is the maximum value of F(y) that occurs in the BL profile and y max is the 
value of y at which F max occurs. In wakes the exponential term is set equal to zero. The function 
Fkleb is the KlebanofF intermittency factor given by 


Fkleb{v) 


1 + 



Cj<JLEBy _ \ 6 
y\iA\ ’ 


(3.29) 


and C'kleb = 0.3. 

Although not seen explicitly, the Reynolds number enters into the computation of eddy 
viscosity through the computation of j/ + . When the variables in equation (3.24) are nondimen- 
sionalized, the following expression is obtained: 


y 


+ 


Rc Pu'T’w - 

— y 


(3.30) 


For more discussion on the Baldwin-Lomax turbulence model, the reader is referred to the original 
paper (ref. 10). 

The distance from the wall y used in the preceding formulations is an easily computed ar- 
clength for points near a simple planar surface. The same lengt h in corner regions was determined 
from a formulation developed by Hung (ref. 65), which took into the consideration the normal 
distance from both walls as follows 
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(3.31) 


d = 


2 y~_ 

[( y + z ) + vV 2 + - 2 ' 


where y and 2 represent the arclengths normal to the intersecting walls. 

I 11 the wake regions the distance from the wall is calculated normal to the zone centerline, 
which hypothetically coincides with the wake centerline (fig. 20). For present purposes the as- 
sumption is probably valid at low to moderate angles of attack up to one half chord aft of the 
trailing edge. Since the shear energy is not large in the wake, compared to the BL, the primary 
improvement in the wake flow calculation comes more from the refinement of the grid, than from 
turbulence model inputs. This use of the wake centerline will be more accurate if the wake grid 
is made adaptive in future code improvements. 

The Baldwin-Lomax model has obvious limitations because of its reliance on the eddy vis- 
cosity concept, but it has proven to be reasonably accurate for predicting moderately separated 
flows at reasonable computational expense. 
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4. STRAKE- WING-FUSELAGE SOLUTION RESULTS 


4.1 Introduction 

The determination of which solutions were to be computed was driven by the aircraft design 
goal: to exploit the formation of the separated vortex from a strake or leading edge extension 
for the generation of so-called “nonlinear lift.” To date, the understanding of these complicated 
interactions has been restricted to observations of experimental wind tunnel and water tunnel 
flow visualization, and evaluation of limited measurements. Wind tunnel visualization has been 
comprised of smoke trails, which dissipate; light sheet illumination, which is planar; and surface 
oil flow traces. Laser velocimetry and anemometry have provided nonintrusive measurements of 
vortical structure much improved over those attainable with probes, but these techniques are 
still too expensive and cumbersome to be used routinely. This is especially true for evaluation of 
dynamic phenomena over extended sampling durations. Although water tunnel flow visualization 
has proven more tractable, its Reynolds number limitation, caused by cavitation effects and dye 
dispersal at higher water-flow velocities, has made questionable the utility of using this technique 
for anything but the most qualitative of comparisons (ref. 66). The computational efforts are 
intended to overcome these deficiencies, but must be compared against experimental results, 
within the limitations of their domain of validity, to verify the computations. 

To this end, cases were solved over a wide range of and Reynolds numbers to get the 
broadest possible comparison with available and attainable experimental data. Solutions at higher 
angles of attack were emphasized to insure the formation of a strake vortex, and when possible a 
wing vortex. Wing tip vortices, which others have fine tuned the TNS code to capture (ref. 67), 
were inadequately resolved in these initial grids. Since the primary interest is the inboard region 
of the flow field, the tip-region resolution is relegated to improvement in future grids. 

Current understanding of three-dimensional separation is principally derived from studies 
about generic shapes such as ogive cylinders. Recent successful computations about compli- 
cated wings (ref. 44) have provided insight into the limitations of the present, solution method in 
properly capturing separated flow. These limitations deal principally with the sensitivity of the 
separation domain to fine tuning of the smoothing parameters, and relaxation of the turbulence 
model in areas of massive separation. There have been some successful computations of aircraft 
configurations, but for relatively simple geometries (ref. 16), with limited separation, and often 
without proper capture of flow physics (ref. 6). 

In this chapter, the computational aspects of the TNS zonal algorithm are presented, and 
then the mathematical and topological considerations of separated flows will be discussed to 
provide a framework for viewing the solution results. The impact of the artificial dissipation 
and turbulence models on separated flow will be considered. Then the experimental test data 
base will be described, followed by presentation of the laminar and turbulent solutions. The 
flow-fields will be evaluated through topological analysis of computed particle traces, comparison 
to wind and water tunnel flow visualization and comparison to wind tunnel pressure distribution 
measurements. 


4.2 Computational Aspects of the Zonal Algorithm 

It will be beneficial to investigate some of the measures of solution efficiency to evaluate 
how the improvements driven by systems considerations affect the duration and expense of the 
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solution process. This study presents the development and application of an engineering tool; but 
the skillful utilization of this tool is essential to getting the most out of it. 

4.2.1 Measures of solution efficiency.- In Chapter 3, the modifications to the block tridi- 
agonal matrix structure which rendered the algorithm scalar pentadiagonal were discussed, along 
with an evaluation of improvements of both efficiency and accuracy because of the easy incorpora- 
tion of a fourth-difference implicit artificial dissipation operator. In this section we will continue 
with an assessment of data management strategies from efficiency, accuracy, and reliability per- 
spectives. 

The convergence rate of the diagonal version of the code versus the block tridiagonal version 
has already been depicted in figure 25 (ref. 50). The geometry used in that sample simulation was 
a NACA 0012 wing at free-stream flow conditions; = 0.826, a = 2° and Re = 8.0xl0 6 . This 
case involved a strong shock with a moderate amount of separation. As such it is representative 
of the efficiency improvements realized in the present transonic solutions of a more complex 
configuration. 

The data management of this large application code, discussed in general terms in Chapter 
2, can now be assessed through evaluation of some sample solution statistics. Memory, speed, 
and I/0(input/output) requirements of the TNS code are quite substantial. The key element for 
handling the memory and I/O requirements is the solid state device (SSD), especially in handling 
the transfer of metrics, as the example will show. 

The metrics, which are stored in three different orientations, constitute a large body of 
information. Since each of the metric arrays is required in main memory several times for each 
grid zone during each iteration, the I/O requirements are formidable. When utilized properly 
the SSD handles this without problem, but if insufficient SSD is requested then the information 
which will not. fit in the available SSD is automatically fed to disk storage. As the computational 
statistics displayed in table 6 demonstrate for a turbulent, calculation of the 18 zone grid (M x = 
0.9, Re = 4.53xl0 6 and a = 4°), the I/O becomes prohibitively expensive both in duration and 
computer cost. To put tangible perspective on the time expense, the disk-storage-intensive run 
took 4 days of intermittent run-time, which needlessly limited the access of others to the system. 
The SSD optimized version ran cleanly in 7 real-time clock hours. 


TABLE 6.- COST STATISTICS FOR TWO TNS PROGRAM SOLUTIONS 
(18 zone wing-body, Moo = 0.90, Re c = 4.53 x 10 s , a = 4°) 


Quantity 

95% storage in SSD 

50% storage in SSD 

Iterations 

500 

500 

CPU Time (Hr) 

2.80 

2.89 

I/O Requests 

2.49xl0 6 

2.51xl0 6 

Memory x CPU time* 

9011 

9321 

Memory x I/O time* 

1158 

93878 

SSD I/O Time (Sec) 

433.7 

66.7 

Disk I/O Time (Hr) 

.20 

9.85 

Cost (computer dollars) 

sk * * * 1 1 * 1 1 

5561 

46919 


* Million words per second. 


The statistics in table 6 paint a sobering picture of the I/O costs that inattention to proper 
employment of SSD storage will cause. The two jobs compared are identical, except that one 


34 



achieved 95% of its temporary storage on SSD because of proper accessing job control language, 
and the other accomplished only 50% SSD temporary storage, the remainder sent by default to 
disk. The 95% solution also had some unnecessary data storage trimmed out, which resulted in the 
slight difference in CPU time, I/O requests and Memory X CPU time. The real difference in the 
two is reflected in the I/O time figures of merit, which show two-orders-of-magnitude difference, 
resulting in a computer accounting cost difference of an order of magnitude. The message clearly 
is that care should be taken to use SSD to the maximum extent possible. 

4.2.2 Vectorization enhancement.- To take advantage of the high-speed characteristics 
of a vector processing systems like the Cray XMP series, it is essential to have an understanding 
of the fundamental abilities of a given Cray compiler as it accomplishes vectorization on the 
innermost do-range of a multidimensional do-loop. Some fundamental considerations include 
(1) that the longest string (i.e., the finest grid dimension) of those ranges available should be 
placed in the innermost loop; (2) dependence on information from a previous increment within 
the innermost loop should be avoided; and (3) “if” statements in the innermost loop should be 
avoided, or if unavoidable should be handled using specially designed ('ray FORTRAN logical 
statements. There are many subtleties to the correct structuring of vectorizable do-loops, but the 
dividends make the effort, worthwhile. A 14% decrease in run time required for all applications of 
the code was realized by enhancing vectorization in just six do-loops of one often-used subroutine. 
This allows the researcher to obtain faster job turnaround. 

It is worth mentioning that it is precisely this vectorization capability which gives the vector 
processors the ability to run at high “million floating point operations per second” (MFLOPS) 
rates. According to the operations guidelines distributed for the CRAY XMP/48 system, well- 
vectorized FORTRAN codes using the standard compiler without any special optimization should 
be able to attain 100 MFLOPS. Since we have clocked 63 MFLOPS for the wing-alone version of 
the TNS code, there seems to be room for improvement. 

4.2.3 Zonal restructuring.- A grid inadequacy which was uncovered during the presenta- 
tion of one of the initial solutions provided the opportunity to demonstrate the versatility of the 
zonal approach for grid modification. This inadequacy was that the grid was too coarse in the 
circumferential direction around the strake edge, resulting in the inability to properly capture 
the vortical effect, and resultant low-pressure region on top of the strake for high angle of attack 
flows. Although the solution converged, the pressure in that region was higher than appropriate. 
The alternatives were to generate a new coarse base grid with finer circumferential spacing, or to 
replace a zoned region of the coarse grid with new zones, which would possess the requisite fine 
spacing; this was the course chosen. The resultant approach replaced the forward fuselage zone, 
which had been designated ZONE 1, with ZONES 1,17, and 18, as depicted in figures 25-26. The 
new ZONE 1 and ZONE 17 represent subsets of the previous ZONE 1, since the coarse circum- 
ferential spacing could be retained in those regions. As will be described later in this chapter, the 
refined spacing in the strake region now covered by ZONE 18 is adequate to capture the most 
sensitive of vortical flow details. This process was accomplished with relatively simple and quick 
modifications to the ZONER grid subdivision program, with results which demonstrate its utility 
for geometry redefinition. 


4.3 Mathematical and Topological Considerations 

4.3.1 Background.- The use of strakes and similar devices which generate nonlinear lift, can 
be characterized as the manipulation of separation effects in a positive fashion. The evaluation of 
the flow fields for two-dimensional applications can be relatively straightforward, since the concept 
of flow reversal, in which flow is opposite to the free-stream direction, provides a boundary. For 
the three-dimensional problem, which is appropriate for aircraft configurations, this concept is no 



longer that helpful since the direction reversed is itself ambiguous. The proper approach to three- 
dimensional flow-field analysis has, therefore, received much attention and resulted in contrasting 
positions on the best way to quantify the resultant phenomena. 

The development of measures of three-dimensional separation have sprung from the physical 
manifestations available from wind tunnel tests. Maskell (ref. 68) was attempting to explain the 
formation of oil streaks on models which had been coated with oil dots prior to tunnel operation 
when he categorized the pattern of skin-friction lines close to the surface as the skeleton structure 
of the entire flow field. 

Legendre (ref. 69) introduced the concept of limiting streamlines which possessed the qualities 
of a continuous vector field. Since streamlines represent the local velocity direction at each point in 
the flow field, they would not exist on a no-slip surface itself, where the velocity is zero. Through 
any nonsingular point there could be only one trajectory. Many trajectories could converge 
to or diverge from the singular or critical points in a variety of fashions. The mathematical 
categorization of these various singular points, and the lines which emanate from or enter them, 
provide the mathematical framework within which to evaluate and quantify three-dimensional 
flow fields. 

It was Lighthill (ref. 70) who tailored this approach to viscous flows by treating a continuous 
vector field of skin friction lines rather than limiting streamlines. This analysis provided validity 
to the oil-streak techniques employed by Maskell (ref. 68), Maltby (ref. 71) and others. Coincident 
with the development of Light hill’s treatment was the shared observation that three-dimensional 
separation was manifested by the convergence of skin-friction lines onto a particular skin-friction 
line. From this line a stream surface leaves the solid surface, and typically coils up in a vortical 
fashion. There is general agreement that this is a necessary condition for separation to occur, but 
its sufficiency is questioned. 

4.3.2 Mathematical foundations for the skin-friction field.- To examine the condi- 
tions for separation and behavior of critical points, it is necessary to consider the wall shear 
stress field. Skin-friction lines are defined as the integral curves of the wall shear stress vector 
exerted by the fluid on the wall. If at any point the wall shear stress vanishes, the field possesses 
a saddle or a nodal point, at which the skin-friction lines’ directions are indeterminate. The 
“phase-plane” and “phase-space” methods of exploring the properties of solutions of ordinary 
differential equations have proved to be extremely successful in the field of nonlinear dynamical 
systems. Oswatitsch (ref. 72), and Lighthill (ref. 70) examined viscous flow patterns close to a 
rigid boundary and classified certain classical points which can occur. They used a mathemat- 
ics which is equivalent to the phase-plane trajectory analysis. Also Perry and Fairlie (ref. 73), 
and Hornung and Perry (ref. 74) applied these techniques to fluid flow problems. The following 
analysis outlines the mathematical procedure for the topography of skin-friction lines and critical 
points. 

The system of equations for particle paths in a steady flow is given by 

= v(r,y,z) 

= v{x,y,z) (4.1) 

= «’(ay V, -) 

where .t, (/, and z are Cartesian coordinates and u,r, and w are the Cartesian velocity components 
along .r, i/, and r, respectively. The skin-friction lines, or limiting streamlines, are obtained by 


dx 

dt 

dt 

dz 

Jt 
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exploring the limiting behavior of equations (4.1). But because of the no-slip boundary condition 
at the wall, equations (4.1) are trivial (i.e., u = v = w = 0, as 2 — > 0; assuming z is normal to the 
surface). At this point, it is helpful to introduce a pseudo-time dr such as 


Rewriting, equations (4.1) yield 


dr — z dt 


x - 


y = 


dx 

dr 

dy 

dr 

dz_ 

dr 


u 

z 

V 


w 


(4.2) 


(4.3) 


It is necessary to analyze the limiting behavior of the shear-stress tensor given in equations (3.3). 
At the surface, the following conditions for the velocity components and their derivatives hold: 


u = v = w = 0 

U x = Uy — V x = Vy = W x = Wy = 0 

but from the continuity equation, w, = 0 also. Hence equations (3.3) reduce to 


(4.4) 


Tex — ^~y y — z 


— r xy — 0 


du. 

dv 

T = »Yz 


(4.5) 


Taking the limits as z — > 0 for the last two of the preceding equations yields 


u T x: 


v 





9(*,y) 


(4.6) 


or 

u r w 

z n 


in which u = («., v) T , t w = (r x: ,T yz ) T , and f(x,y) and g{x,y ) are nonlinear functions. 


(4.7) 


It is assumed that, for steady flow of an incompressible Newtonian fluid with constant vis- 
cosity at. finite Reynolds number over a smooth continuous wall, the velocity and pressure are 
regular, so that local solutions can be written as Taylor series expansions in the space coordinates. 
The mathematical support for this assumption is discussed by Ladyzhenskaya (ref. 75). Equa- 
tion (4.7) reveals that the vector quantity 11 / r is equal to r u ,//< , and the vector field 11 fz lias the 
same integral curves as the wall shear stress, i.e.. the skin-friction lines. When u / z is expanded 
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into a Taylor series about the critical point, and only the lowest order terms are retained, the 
following linear system of first-order ordinary differential equations is obtained: 


X 


a b 


X 

y . 


c d 


y . 


where a.b,candd are constants of the expansion, or 

x = F • x 


(4.8) 


(4.9) 


Hence, the assumption that the velocity components are Taylor series expandable about a critical 
point is equivalent to assuming that the (autonomous) partial differential equations for fluid flow 
can be reduced to a set of ordinary differential equations (ref. 73). The matrix F of equation (4.9) 
has eigenvalues Aj and A 2 which may in general be real or complex. The corresponding eigenvector 
slopes are 

mi = (Aj — a)/b - c /{ Aj - d) 

(4.10) 

tt?2 = ( A > — a)/b — c/( A 2 — d) 

which, for the case of A] and A 2 real, may be shown to correspond to the slopes of certain 
trajectories, the sepratrices, which emanate from the critical points. 

There is just one skin-friction line through each point on a surface, except a point of separat ion 
or attachment, where r w = 0. These points are the singular points of the differential equations 
governing the topography of the skin-friction lines. Such singular points are classified (ref. 76) 
depending on the values of the following quantities: 


p = —(a + d) = —trace F 

(4.11) 

q = (ad - be) = det F 

and 

Ai, 2 = ~[pT(r-4g) 1/2 ] (4.12) 

The classification of possible critical points is shown in figure 29 from reference 73. A singular 
point where q < 0 is a “saddle point.;” a point where q > 0, however, is a “nodal point.” A nodal 
point (fig. 29(a)) is a point where all of the skin-friction lines except one (labeled A A in fig. 29) 
are tangential to a single line BB. A focus (fig. 29(b)) differs from a nodal point in that it has no 
common tangent line. An infinite number of lines spiral around the singular point, either away 
from it (a focus of attachment.) or into it (a focus of separation). At a saddle point (fig. 29(c)), 
there are only two particular lines (C'C and DD) that pass through the singular point. These lines 
represent eigenvectors. The flow directions on opposing sides of the singular point are inward on 
one particular line and outward on the other line. All of the other lines pass by the singular 
point in directions consistent with the directions of the adjacent eigenvectors. Critical points 
corresponding to values of p and q along the axes (p — 0 and/or q — 0) and on the parabola 
p = 4q 2 are “degenerate” forms. 

4.3.3 Topography of skin frict ion lines.- Singular points have certain characteristics 
that largely determine the distribution of skin-friction lines on the surface. The nodal point 
of attachment is typically a stagnation point on the forward face of the body, where the free 
stream attaches itself to the surface. Hence the nodal point of attachment behaves like a source 
from which the skin-friction lines issue and circumscribe the body. The nodal point of separation, 
however, behaves like a sink where the skin-friction lines on the body surface terminate. I 11 other 
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words, this is the rear stagnation point of the body. The saddle point acts typically to separate 
the skin-friction lines issuing from these nodes. The total number of singular points for a possible 
pattern on a smooth surface is subject to a topological rule that the number of nodal points must 
exceed the number of saddle points by 2 (ref. 76). A physical explanation for this is presented by 
Lighthill (ref. 70): since there are infinite number of skin-friction lines on the surface, and they 
must begin and end somewhere, there is at least one nodal point of attachment and one nodal 
point of separation. If there are two nodal points of attachment, however, the skin-friction lines 
from each node must somewhere run into each other; this necessitates the introduction of a saddle 
point in between. 

This type of an arrangement, which has been detected on the forward facing canopy of the 
Space Shuttle, is depicted in figure 30 taken from (ref. 69). This same combination of two nodal 
points of attachment and one saddle point of separation is quite possibly present on the F-16A 
canopy at low angle of attack, but such a feature will not be resolved until the grid is made finer 
in that local flow region. 

If there are n nodal points of attachment, there will be (n — 1) saddle points accordingly. 
Hence the whole combination is equivalent to a single node of attachment which behaves like 
a “source” as explained before. Similarly, m nodal points of separation and (m — 1) saddle 
points behave like a “sink” into which the skin-friction lines emanating from the “source” vanish. 
Therefore, there are (m + n) nodal points and (m + n — 2) saddle points, and the topological law is 
satisfied. A particular skin- friction line emerging from a saddle point prevents the lines from the 
nodal points of attachment from crossing each other. This particular skin-friction line is called 
a line of separation. Skin-friction lines from either side asymptotically converge on the line of 
separation. Hence, it is possible to define the line of separation as a particular skin-friction line 
toward which other skin-friction lines converge rapidly (in an asymptotic manner, but not a cusp- 
like manner (refs. 70 and 73). Conversely, a line of reattachment. is defined as a particular line 
from which other skin-friction lines diverge rapidly. It should be kept in mind that this definition 
is associated with the necessary condition of the separation. Whether it is also sufficient continues 
as a matter of debate (ref. 77). 

In the vicinity of the line of separation, the limiting streamline must leave the surface. This 
can be shown as follows: at a very small distance z from the surface the velocity was found 


z n 

Continuity must be satisfied in a stream tube, for which the volume flow rate M is expressed as 

M — z d ^ = const. (4.13) 

From equation (4.7) 

q=-zr u . (4.14) 

where d is the distance between the streamlines, q = ] u j and r !( . — j r„, j. From equations (4.13) 
and (4.14) 

- - (dr„.) _1/ ' (4.15) 

Hence, the height of the limiting streamline c above the surface grows rapidly as the line of 
separation is approached. According to equation (4.15). the mechanism which causes the limiting 
streamlines to leave the surface is twofold: first, the resultant skin-friction t w becomes vanishingly 
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small near the saddle point, and second, the distance d between the adjacent limiting streamlines, 
falls rapidly as the limiting streamlines converge toward the line of separation. 


4.4 Effects of Artificial Dissipation and Turbulence Models 

Turbulence and artificial dissipation models play crucial roles in the simulation of separated 
flows and shock/BL interactions. Since there is an intimate relation between the skin friction field 
and flow separation, correct computation of BLs is very important. To assess the relative effects 
of these models on the BL characteristics, a series of test cases was computed by Kaynak et al. 
(ref. 44) for attached and shock-induced separated three-dimensional flows. The extrapolation of 
the conclusions drawn about the utility of dissipation and turbulence modeling concepts for the 
strake- wing-body computed herein is straightforward. 

Kaynak compared computed results of a two-dimensional att ached flow around a NACA 0012 
airfoil at = 0.5, Re = 2.89 xlO 6 and a = 0° to experimental wind tunnel test data taken 
by Thibert et al. (ref. 78) for pressure coefficient distributions, and found good agreement. He 
assessed BL capture abilities of the TNS code by comparison against computed profiles from a 
proven transonic airfoil, interactive BL code written by Steger and Van Dalsem (ref. 79) called 
TRIVIA. It was determined from the comparison that the TNS code was not developing correct 
BL profiles. The smoothing scheme in use with the diagonalized algorithm was developed for an 
inviscid solver (ref. 57). With this scheme the smoothing varies inversely with grid cell-volume, 
which is inappropriate for a viscous grid. Kaynak corrected this problem by specifying a Mach 
function: 

{ A/ 2 , for attached flows 

(3.32) 

M, for separated flows 

to be multiplied by the original smoothing operator indicated in equation (3.19). This approach 
reduces smoothing for both separated and attached subsonic flows, but differentiates between the 
two. The amount of correction was determined computationally by trial-and-error, and in all 
cases forces the artificial viscosity to zero at solid no-slip boundaries. This correction has not 
been employed in the present solutions, but should be applied when accurate coefficient of drag 
(C’a) calculation is desired, or for BL profile portrayal. 

Kaynak et al. (ref. 44) also found that modifications to the turbulence model were necessary 
to properly capture separation effects. The Baldwin-Lomax model, which is an equilibrium model, 
predicts too sharp an increase in the outer layer eddy viscosity through the separation shock. The 
increase of eddy viscosity from its local equilibrium value far from the wall leads to more turbulent 
mixing and to more shear st ress to balance the adverse pressure gradient , suppressing the tendency 
toward separation. Rotta (ref. 80) concluded from experimental data that when turbulent flow 
is perturbed from its local equilibrium state, a dist ance of about one order of magnitude greater 
than the boundary layer thickness is required to attain a new equilibrium state. To account 
for the upstream turbulence history effects, Shang and Hankey (ref. 81) used a relaxation eddy 
viscosity model; i.e., 

( = (cq~ (f cq ~ (o k*P ( ~ ' r - A ' r ~ ) (3.33) 

where e is the turbulent dynamic eddy viscosity, t cq is the local equilibrium eddy viscosity 
evaluated from equation (3.20), Co is the eddy viscosity at the upstream location .r 0 , and \ is 
the relaxation lengt h. A good review of turbulence model relaxation techniques can be found in 
Hung (ref. 82). Conceptually, equation (3.33) approximates the experimental observation that, 
in an abrupt disturbance of a turbulent flow, the Reynolds stress remains nearly frozen at its 
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initial value while it is being convected along streamlines, and then exponentially approaches 
a new equilibrium state. In a numerical calculation, the location from which the relaxation 
process is initiated, x 05 and a relaxation length scale which describes the exponential decay of 
the eddy viscosity distribution, A, must be specified. There are two limiting cases which bound 
the relaxation length. For A = 0, the turbulent eddy viscosity equals the local equilibrium value, 
and for A = oo, the initial value e 0 is frozen and is used everywhere in the region x > x 0 . Via 
numerical experimentation, Kaynak found A = 40 to be a good value. 

His implementation of the relaxation turbulence model was applied to the computation of 
three-dimensional separated flows on wings for the first time, and improvements were observed 
in the skin-friction topography as well as the pressure coefficient. This approach has not yet 
been applied to the solutions of F-16A flow-fields, in part because the supercritical wing design 
delays massive separation to high angle of attack. When these high angle of attack solutions are 
obtained, this technique and its impact on the separated flow field should be evaluated. 

A special consideration with the application of the Baldwin-Lomax (ref. 10) turbulence model 
for a wing-fuselage is the proper capture of cross-flow separation effects on the leeward side of 
the fuselage at high angle of attack. Cross-flow separation, as depicted in figure 31 taken from 
Degani and Schiff (ref. 83), occurs when the flow from windward to leeward separates into a coiled 
vortex. This primary vortex reattaches near the symmetry plane, often generating a secondary 
vortex inboard. The proper capture of this secondary vortex, and even possible tertiary vortices, 
is essential to correctly reflecting the physics driving the primary vortex, in its formation and 
shedding process. The strength and impact of these vortical constructions increases with angle 
of attack. 

The proper capture of these constructions with the Baldwin-Lomax turbulence model is 
problematical because of the manner in which the function F mai (See Sec. 3.6.3) and subsequently 
(Pt)ouUt (eq. 3.25) are determined. The difficulty stems from the calculation of F(y) (eq. 3.28), 
for which the maximum is determined on a ray which extends from the surface through the BL, 
as shown in figure 31. For a simple BL profile a well defined peak will be easily determined, 
for instance at location b in figure 31. For cross-flow separations, two or more peaks can occur 
because of the distinct vortical constructions (fig. 32), the inner one (designated a in fig. 31) of 
which gives the proper F(y) value. If the maximum F(y) for the higher peak of the outer primary 
vortex is incorrectly utilized, the secondary and tertiary aspects of the flow field may be washed 
out. For the F-16A geometry of the present study, secondary and tertiary constructions were 
captured without including the Degani-Schiff F(y) correction. The reason for this unexpected 
success is unclear, except that the sharp edge of the strake (fig. 27) may cause this geometry to 
be less susceptible to the problem. It is likely that in fuselage regions upstream of the strake this 
problem will be significant at higher angle of attack, requiring the correction outlined in (ref. 83). 


4.5 Simulation of F-16A Faired-Inlet Flow-Fields 

Both laminar and turbulent calculations were performed because there was a variety of 
experimental test media available to verify the computational results. Turbulent flow-fields were 
calculated at = 0.5, a = 0° and at = 0.9 for 2°, 4°, and 6° angle of attack at. a Reynolds 
Number (Re) of 4.5 xIO 6 (based on wing-root chord), to match the wind tunnel data. The laminar 
solutions were computed at Moo = 0-3 and 0.6, with 10° angle of attack and Re of 1.5 x 10 4 , 
to match the waler tunnel Reynolds number, and to evaluate a possible dependency of the 
strake vortex calculation. A quasi-laminar case was also run at the wind tunnel Reynolds number 
to determine Reynolds number effect, on vortex capture. This case also attained a higher angle 
of attack than we could initially achieve with the finer grid required for turbulence modeling. 
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4.5.1 F-16A test data.- The primary experimental data used to confirm the accuracy of 
the computed results was from a 1975 wind tunnel test directed by Hammond (ref. 39) using a 
one-ninth-scale F-16A model. This test was accomplished in the USAF AEDC 16T at a standard 
Reynolds number of 2.5 x 10 6 per foot, covering the range from 0.60 to 1.55. The angle of 
attack ranged from -5° to 28°, and for this model, tunnel blockage is estimated to be less than 
one percent. Force balance load measurements were taken, as were pressure measurements from 
510 orifices located on the model. Transition grit was affixed to the model in a shallow ring 1.67 
inches aft of the nose, and within one-tenth inch of the leading edge of the wing and tail surfaces 
to promote correct turbulent BL generation. Grit was also placed 0.73 in. aft. of the inlet lip on 
the fuselage underside for the same purpose. 

Oil flow visualization was available from a 1972 YF-16A test entry, which used a geometry 
equivalent to that of the 1976 F-16A test. The visualization was accomplished by painting the 
model with oil and then varying the angle of attack at constant M^. Apparently the angle 
of attack was maintained briefly at a given setting to stabilize the oil flow pattern prior to 
photographing. The angle of attack was then increased while the tunnel continued to run, without 
applying more oil. This “alpha sweep” technique may result in hysteresis from a lower angle of 
attack oil streak pattern affecting the position and thus shapes of subsequent patterns at higher 
alpha. There is also the possibility of the oil setting up and resisting change. For these reasons 
oil flow visualization must be considered as an aide to the understanding process. It should not 
be relied upon to provide a definitive characterization of the surface streamline behavior. Oil flow 
visualization was also available from a 1972 Pre-YF-16A geometry evaluated in a wind tunnel 
test run at the Cornell Aeronautical Laboratory facilities, Buffalo N. Y. Since this configuration 
possessed an unreflexed strake (fig. 8), comparison of strake-related flow phenomena would be 
misleading. However, general flow features exhibited by oil flow on the wing are still of interest. 

Water tunnel tests were conducted on a one-seventy-second-scale F-16A plastic model in 
the NASA Ames Dryden Research Center Flow Visualization Facility (FVF). The test section 
in that facility measures 16 in. by 24 in., and is oriented vertically, with plexiglas sidewalls for 
all aspect viewing around the longitudinal centerline. Typical Reynolds number ranges varied 
from 1.5 x 10 4 to 2.5 x 10 4 per foot, with typical freestream fluid velocity of about 0.25 ft /sec for 
optimum vegetable dye-trace flow visualization (ref. 84). Both extensive still photography and a 
VCR videotape were made of selected visualization features for extended viewing of significant 
aspects of the flow. 

Evaluation of the computed data is greatly facilitated by use of graphics workstations and 
post processing graphics software entitled “PLOT3D” developed by Peter Bulling of Ames Re- 
search Center (ref. 85). A variety of scalar and vector functions may be presented in three- 
dimensional visual perspectives, as well as computed particle traces and special surfaces of inter- 
est. 


As discussed by Kaynak et al. (ref. 44) the particle traces can be divided into two groups: 
“st reamlines” and “sect ional streamlines.” These two terms are explained in figure 33. If a particle 
is released from point A, its path may be plotted in several ways, including the three shown, 
AB, AC, and AD. If its motion is not confined to any surface, it moves according to the fluid 
dynamic forces and becomes a “streamline.” Note that in steady flow pathlines are equivalent to 
streamlines. If the motion of the particle is confined to r - y or y — r planes, “sectional streamlines” 
are obtained; the lines AD and AC, respectively. It is not hard to imagine a sectional streamline. 
For example, the line AD may be likened to the motion of a free-surface fluid particle when a 
body is immersed into a liquid and set into motion. In this case, the motion of the particle is 
confined to the free-surface of the liquid. The sectional streamlines are computed from only the 
components of velocity on that particular surface; i.e., the projections of the velocity vectors on 
the surface. Note that the sectional streamlines are not the same as the projection of streamlines 
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on the surfaces. For example, if the streamline AB was projected on the x — y or y — c surfaces, 
lines other than AC or AD would be obtained. 


4.5.2 Turbulent flow cases.- Solutions have been run at = 0.9 for angles of attack 
of 2°, 4°, and 6° at a Reynolds number of 4.53 xlO 6 , with the convergence criterion being at 
least a three-orders-of-magnitude drop in the L2 Norm of the residual in all zones. Typically 
the drop ranged from three to five orders of magnitude from the viscous to the inviscid zones. 
Of these solutions the most difficult to compute was the higher angle of attack case: at lower 
alpha the separation at the strake edge and the gradients at the wing leading edge are not as 
pronounced; therefore, the time step could be increased and a converged solution was achieved 
more rapidly. The Moo = 0.9, a = 2°, 4° and 6° cases took approximately 2000, 3000, and 
3200 iterations, respectively, to converge. To give perspective, the 3200 iteration requirement of 
the — 0.9, q = 6° case cost 16 hours of Cray XMP central processing unit (CPU) time. 
Since efficient programing has not been carefully attended to in these proof-of-concept solutions, 
it is expected that the process can be speeded up considerably with fundamental corrections to 
enhance vectorization and minimize time consuming data transfer. Preliminary indications point 
to a near term savings of 50%. 

A 6° angle of attack is not high relative to what the F-16A can achieve in flight, but is high 
enough to generate the fluid phenomena of interest, and is the maximum alpha we have been 
able to compute for these turbulent flows to date. Grid clustering in the directions normal to the 
aircraft surface provide adequate definition of the BL. The first grid point off the surface in the 
viscous sublayer is typically at a y + value of 5-8 as computed using equation 3.24. 

As shown in figure 34, the computed results are quantitatively in good agreement with the 
experimental pressure coefficient distributions on the fuselage centerline and at four spanwise 
stations on the wing. The compression and expansion over the canopy, the suction peak near the 
leading edge, and the oblique/normal shock combination near the wingtip are all captured. Finer 
streamwise resolution, to be used in forthcoming grids, is needed to improve the accuracy of their 
depiction. Similar trends were observed in the comparison of the Moo = 0-9, o = 4° computed 
data to experimental data, as depicted in figure 35. Preliminary computations (ref. 22) indicate 
that with increased resolution, the computed aft normal shock will extend farther inboard, in 
better agreement with experiment. The shock features of this solution are well detailed by the 
computed surface pressure contours and Mach contours near the BL edge shown in figures 36 and 
37. It is worth mention that the Mach contour clustering indicated at the wing-fuselage juncture 
region in figure 37 is a result of graphical presentation technique. It is caused bv the use of 
constant-index normal surfaces, which dip into the BL in corner regions. 

Experimental oil-flow visualization for these conditions (Ref. 20) is shown in figure 38, while 
figure 39 shows a enlargement of the strake- wing juncture region. Comparisons of figures 38 and 
39 with the analogous computed surface-particle traces in figures 40 and 41 reveal agreement that 
the flow field, on the wing upper surface, remains attached from the leading edge to the shock- 
induced separation. On the strake the spacing in the spanwise direction, approaching the edge, 
has been refined with additional zones in subsequent grids. The computed secondary separation 
and reattachment lines are well captured, and compare favorably to the oil flow. The very faint 
separation line near the inboard leading edge in figure 39 is not captured in the computation; 
inadequate leading-edge definition may have caused the discrepancy. The experimental oil flow 
(fig. 38) also shows the normal shock located farther aft than does the computed result, which 
is appropriate since wind-tunnel walls were not modeled. As noted with the pressure coefficient 
distributions, the computed normal shock should also extend farther inboard. Comparison be- 
tween the experimental and computed results in the vicinity of the wingtip is difficult, since the 
wind-tunnel model had tip missiles and missile rails, whereas the computed configuration does not 
yet contain these features. The spanwise resolution is also inadequate to accurately capture tip 
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separation until more points are added, as was done by Srinivassan et al. (ref. 23) in a tip- vortex 
study. 

Figures 42 and 43 display contours of constant vorticity and total pressure in the cross-flow 
plane located just aft of the strake-wing juncture. The perspective is from upstream, slightly 
inboard and above. In both representations the core of the primary vortex is well defined, as is 
the shear energy in the BL near the wing leading edge. The corresponding cross-flow velocity 
vectors are shown in figure 44. Examination of the core region reveals the clockwise flow direction 
of the vortex. The discontinuities in the vorticity contours (fig. 42) occur at a zonal interface, 
and are caused by one-sided spatial differencing used to calculate this derived quantity from 
the computed velocities. For comparison, no discontinuities are apparent in the total pressure 
contours (fig. 43), since they are obtained directly from algebraic manipulation of the primitive 
solution variables. For these conditions the secondary vortex suggested by the surface-particle 
traces and oil flow patterns is compressed well within the BL. 

These results may be compared to those of Karman et al. (ref. 6), in which a very accurate 
surface modeling of the F-16A geometry was fit with a coarse Euler grid of 500,000 points. 
Unfortunately the explicit Euler approach employed was extremely expensive, costing over 20 
hours of CRAY XMP/48 CPU time without converging an = 0.9, a = 4° case. 

4.5.3 Laminar flow cases.- A quasi-laminar subsonic solution was obtained for = 
0.6, a = 10°, and Re = 8.0 xlO 6 . Although the Reynolds number is well above the expected 
transition Reynolds number, the turbulence model was not employed because it would not have 
been effective with the slightly coarsened grid required to converge the solution in this higher 
incidence case. It took 1200 iterations to drop the error residual over three orders of magnitude 
in all zones, at a cost of 6 cpu hr of Cray XMP/48 time. Subsequently, laminar cases were 
computed using Re = 1.5 xlO 4 to match water-tunnel test conditions for comparison with flow 
visualization (Ref. 21). 

Cross-flow velocity vectors for this case are shown in figure 45. The viewpoint is the same 
as that used in figure 44. The strake-generated vortex appears to be located farther from the 
strake-wing juncture surface than was the case for Moo = 0*9, a = 6° (fig. 44). This vortex also 
has a more visible interaction with the high-velocity flow over the wing leading edge, and with the 
BL. Close examination of the velocity vectors underneath the clockwise-rotating primary vortex 
reveals a counter-clockwise-rotating secondary vortex, and adjacent to it a clockwise-rotating 
tertiary construction. A diagram depicting the relationship between these vortical structures is 
shown in figure 46, which was adapted from a drawing by Luckring (ref. 24) for an experimental 
strake-wing investigation. The tertiary vortex under the secondary was not detected in the 
computed results, but evidence of its existence in the experimental flow is suggested by the gap 
between the secondary separation lines indicated by the oil ridges observed in figue 39. The 
orange-magenta-colored velocity vectors in figure 45 suggest the interaction of the high-velocity 
flow over the leading edge which is further accelerated by the vortical strake-flow. 

Contours of vorticity and total pressure for the same solution and perspective are displayed 
in figures 47 and 48. respectively, and yield more information and corroboration of the vortical 
activity present. The blue vorticity contours in figure 47 indicate the presence of the secondary 
vortex core. They change sign from the primary contours above them, highlighting their opposite 
sense of rotation. Although the total pressure contours of figure 48 do not indicate the secondary 
vortex core, the primary core stands out, as does the genesis for the (clockwise-rotating) wing-root 
vortex adjacent to it, indicated by the closed magenta contours. 

Both Wortman (ref. 25) and Erickson (ref. 26) have warned against relying on water-tunnel 
flow visualization for comparison to high Reynolds number experiments and solutions. Erickson 
in particular suggested that low-angle-of-attack comparisons were invalid because the entrainment 
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of turbulent fluid from the wake causes premature vortex bursting. He suggests that for higher 
angles of attack, say above a = 10°, and for higher leading-edge sweep angles the comparisons 
have more validity. In this regime the flows are vortex-dominated, and BL effects, which are 
Reynolds-number- sensitive, are not as significant. For these reasons the following comparisons 
have been used as a guide to understanding the interaction process, rather than as a definitive 
indicator of the flow behavior at higher Reynolds numbers. 

Computed particle traces for this quasi-laminar case are shown in figure 49. The particles are 
released at. selected points near the body surface, and are allowed to move freely in all directions. 
Figure 49 shows that, the strake vortex, indicated by yellow traces at the core surrounded by red 
traces on the periphery, is pulled spanwise underneath the wing- root or wing leading-edge vortex, 
which is represented by a blue core circled by magenta traces. This representation compares 
favorably with the water-tunnel flow visualization shown in figure 50, which was obtained by 
Malcolm and Skow (Ref. 27) for an F-16-like geometry with an extended fuselage and a more 
highly swept wing. 

A schematic examination of vortical interaction in the cross-flow plane, as drawn in figure 51 
for two stages in the evolution, suggests why the strake vortex A is forced under the wing- 
root. vortex B, and why both move outboard spanwise. This perspective is from downstream 
and above the trailing edge, and represents what, might develop from the initial interaction of the 
vortices near the leading edge, to where they might wind around each other near the trailing edge: 
Strake vortex A feels the downward force of wing-root vortex B, combined with the spanwise and 
downward force of fuselage mirror-vortex A’. Vortex B feels the upward force of vortex A, along 
with the spanwise force of mirror-vortex B’. Other mirror effects are deemed to be secondary 
because of distance from A or B. It is also probable that the strake vortex is weakening because 
of the distance traveled from its formation and resultant dissipation, whereas the wing-root vortex 
is both nearer to its source, and is fed by a stronger source. 

To gain further insight, into the behavior and interaction of the strake and wing vortices, a 
water-tunnel test was carried out in cooperation with G. Malcolm and L. Schiff (ref. 21) in the 
Ames-Dryden Research Center’s FVF. The model used was an F-16A, with the exception that 
the fuselage was cut off ahead of the empennage area (figs. 52 and 53). The water-tunnel model 
had a flow-through inlet capability, and tests were conducted over a range of inlet-flow conditions. 
The experimental dye traces more closely matched the computed particle traces when the inlet 
flow was completely blocked, better simulating the pressure build up caused by the faired-over 
inlet.. The comparison to computation is still not that, good, suggesting that, the inlet and aft. 
fuselage differences are significant- Figure 53, which shows the flow-through inlet, result, gives no 
indication of vortex interaction. 

A laminar solution was obtained at = 0.3, o - 10°, and Re = 1.5 x 1 0 4 , to more 
closely match the water-tunnel test conditions. Particles were released from the same locations 
as those for the quasi-laminar case shown in figure 49, and the resulting particle traces are shown 
in figure 54. A comparison of figures 49 and 54 shows that in the lower Mach/Reynolds number 
computation the the traces stayed farther inboard as they moved aft. The wing-root vortex 
still drew the strake vortex beneath it, but for this lower Reynolds number solution they coiled 
prior to reaching the trailing edge, in contrast, to the previous case. An additional solution was 
computed at the lower Reynolds number, Re = 1.5 x 10 4 , but with M^ = 0.6. Particles were 
again released from the same locations as those in figures 49 and 54, and the resulting particle 
traces are shown in figure 55. The Mach number increase appears to shift the particle traces 
farther inboard spanwise, but does not. modify the basic interaction between the strake-generated 
and wing-root, vortices. 
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Other water-flow visualization studies display a significant vortex interaction for similar 
configurations. One is a water-sled hydrogen-bubble flow visualization obtained by Thompson 
(ref. 28), shown in figure 56, for a double-delta wing at slightly higher angle of attack. The 80° 
sweep angle of the first delta closely approximates the the F-16 strake varying-sweep angles, and 
the 40° sweep angle of the second delta matches the F-16 leading-edge sweep. In this flow the 
cores do not translate as far outboard spanwise because there is no fuselage to aid the movement. 
Flow visualization for an F-16A geometry at a = 8°, obtained by Erickson (ref. 26) and shown 
in figure 57, suggests entrainment of the strake vortex by the wing-root vortex. This flow visual- 
ization may be misleading, however, because it appears that the dye may have been released on 
the periphery of the strake vortex core. 

Although the water-tunnel flow visualizations described in the aforementioned experiments 
do not rigorously match the behavior of the computed particle traces, the same fundamental 
vortical convection and interaction processes are captured by both, with variations because of 
Reynolds number, Mach number, and configuration differences. The differences in the vortex 
interactions of the flow-through versus the blocked-inlet experiment (figs. 52 and 53) highlight 
the importance of incorporating the flow-through inlet, with proper mass flow conditions. 

Surface-particle traces are the computational analogy of experimental oil-flow visualizations, 
since the particles are released just above the body surface and move with the local flow velocity. 
Surface-particle traces for the solution computed at = 0.6, a = 10° and Re = 8.0 xlO 6 
are displayed in figure 58, while figure 59 shows the associated critical-point analysis. Figure 60 
depicts an enlarged detail of the strake- wing juncture area of the computed particle traces shown 
in figure 58. It is noteworthy that the separation line near and parallel to the wing leading edge, 
detected in the oil-flow visualization of figure 39, is captured in these computed particle traces. 
Figure 61 shows a topological analysis of figure 60, emphasizing the appropriate secondary vortex 
separation and reattachment line relationships. 

The cross flow separation/Teattachment lines of figure 60 seem to disappear just aft of the 
strake-wing juncture; this suggests the area in which the vortices might physically lift off the 
surface, though they maintain proximity to the wing until they dissipate further aft. The relative 
location of these separation and reattachment lines is the same as displayed for the turbulent 
calculation at = 0.9, a = 6° in figure 39, but in that case the phenomena are pressed closer 
to the strake edge. Whether these cross-flow separations are “true” separations, or something else 
since they neither begin nor seem to end at specific critical points is a hotly debated topic. How- 
ever, their existence in both computational and experimental observations, and the consistency 
of their physical interactions with the surrounding flow-fields is not questioned. 

The relationship of the two counter-rotating foci evident in figure 59 is typical for a low r - 
aspect-ratio wing at moderate angles of attack. This wing analysis is very similar to a delta wing 
analysis done by Legendre (ref. 90) with the addition of a node at the tip reflecting the wing-tip 
vortex origin not present in the delta. The oil visualization photograph of figure 62 for the Pre- 
YF-16A configuration at higher Moo and the same angle of attack shows topological similarity to 
the analysis in figure 59. The same two counter-rotating foci appear as in figure 59, but further 
aft near the trailing-edge flap. This pair produces what has been termed a “musliroom-cap”-type 
separation, as suggested by the free-stream flow which is forced over the separation region which 
they anchor. 

The schematic in figure 63 displays the fundamental fluid interactions occuring above the 
wing fuselage for this solution. The computational analysis at the wingtip is largely estimated 
since the grid is not yet fine enough to resolve the flow detail in this area. A videot ape of the dye 
movement, made during the water-tunnel test, has proven invaluable in developing perspective 
for what has been computed. 
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It is striking that the surface- particle traces for the lower Reynolds number cases; (Re 1.5 x 
10 4 ), Moo = 0.6 (figs. 64 and 65) and Moo = 0.3 (fig. 66) do not display the secondary cross-flow 
separation or reattachment lines on the strake seen in figure 60. The surface-particle traces for 
these cases all flow spanwise to the strake edge to support the primary vortex formation. The 
other wing flow critical point features seem similar to the high Reynolds number result . 

As suggested by Peake and Tobak (ref. 91), a three-dimensional figure which has had critical 
point analysis completed over its entire surface should follow certain rules of topological grammar. 
For instance the total number of nodal singular points (including foci) on a body should exceed the 
total number of saddle singular points by exactly two. This sort of analysis should be accomplished 
on the present configuration once the grid refinement is adequate to capture the major features 
expected. This type of evaluation is helpful in establishing the accuracy of a computed solution, 
in terms of the probability of its correctly capturing the major flow features present. Further 
insight into critical point analysis of flow visualization phenomena is contained in other papers 
by Tobak and Peake (refs. 92-94). 
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5. CONCLUSIONS 


5.1 Summary 

The accomplishments of this research entail advancement in three areas of inquiry: (1) zonal 
solution approaches, (2) application of numerical methods for solving the Euler /Navier-Stokes 
equations, and (3) fluid mechanics evaluation. 

Zonal approaches designed to tailor solution areas for both geometry simplification and the 
application of alternative equation sets selectively have been used before. However, never to the 
author’s knowledge, has this been accomplished for highly complicated fighter-type geometries 
with any form of the Navier-Stokes equations. The thin-layer Navier-Stokes equations were solved 
for a flow-field grid around the strake- wing-body configuration modified from a General Dynamics 
F-16A aircraft geometry. This grid contained approximately 300,000 points. The difficulty of 
stretching grids to properly capture turbulent properties in the BL without needlessly refining 
the grid away from the body has been the major obstacle. The versatility of this method was 
demonstrated by the addition of extra zones to the strake region to improve vortex capture 
during the course of the research. The ZONER approach accomplishes this task in a relatively 
straightforward way which holds promise for application to virtually any vehicle geometry. 

In these applications viscous effects were computed normal to all aircraft surfaces, which 
meant viscous solutions in two directions in corner regions such as the wing-fuselage junctures. 
Viscous effects were also computed in wake zones, including regions behind the wing trailing edge, 
and regions outboard of the tip where the roll-up vortex may be captured. Laminar solutions 
were generated for qualitative comparison with water tunnel flow visualization, and for critical 
point analysis. Turbulent solutions were generated at moderate angles of attack for quantitative 
comparison to wind tunnel data, and for qualitative comparison to wind tunnel flow visualiza- 
tion. The turbulent calculations employed an algebraic eddy-viscosity turbulence model, which 
provided reasonable separation effects when compared to wind tunnel data. 

Extensive use of video graphic simulation was made to analyse the particle trace behavior of 
several of the solutions. This was done to distinguish the primary fluid mechanisms interacting 
in the flow field, to better understand the process in the hopes of using design to control those 
interactions. Various grammars of topological analysis were considered in developing the critical 
point evaluation which resulted. The author has confidence that fundamental interactions which 
occur between the strake-generated vortices and the wing-generated vortices at moderate angles 
of attack have been correctly captured and depicted. 

The proper exploitation of the capture and analysis of separated viscous phenomena was 
the true objective of this research. Good quantitative and qualitative agreement was observed in 
the comparison of the accomplished computed data to both wind tunnel and water tunnel data 
and flow visualization, such that confidence can be placed in the TNS solution approach. This 
tool may now be fruitfully employed to design strake- wing-fuselage configurations for full-aircraft 
evaluations. In modern high-performance aircraft designs, the capture of transonic separated 
aerodynamics is probably the most challenging aspect. In this domain the TNS approach is 
both necessary and adequate. At lower angles of attack, and airspeeds well away from the sonic 
condition, less complicated design tools may be employed. 
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5.2 Recommendations for Future Work 


There are a range of recommendations which may be made for this research covering the 
three areas of accomplishment mentioned above: (1) zonal solution approaches, (2) application 
of numerical methods for the Euler/Navier-Stokes Equations, and (3) fluid mechanics evaluation. 

There is a need to add more points to the grids presently in use to resolve regions of separated 
flow in a lateral sense. The recent availability of the CRAY 2 computer system will allow for larger 
grids without the complicated data transfer problems which previously have accompanied such 
an increase. In addition to resolving the flow features which occur in the present grid, there is 
a strong desire to model the entire F-16A geometry by gridding the tail surfaces and the flow- 
through inlet and exhaust. This will provide more reasonable comparison to the flow-through 
inlet wind tunnel data which is available. It will also make the TNS code a more appropriate 
design tool since propulsion-integration effects are so significant to the inboard aerodynamics near 
the strake-wing juncture. Other researchers are currently at work integrating these features. 

Modifications to the numerical algorithm, other than improvements to its efficiency which are 
evolving during the course of other work, might include the exercise of time accurate capability, 
better capture of shock properties and improved turbulence modeling. The use of the time 
accurate feature of the Beam- Warming (ref. 28) algorithm currently contained within the TNS 
code is problematical because without the use of local time-stepping which destroys time accuracy, 
the code is virtually unusable because it becomes too slow. Even for the steady-state solutions 
presently computed, an increase in efficienty of several factors is needed to make it a practical 
design tool. The incorporation of total- variation-diminishing (TVD) type schemes to improve 
shock capture has been developed for the ARC3D (ref. 27) solution kernel of the TNS Code, and 
may be relatively easy to incorporate. 

The application of the Degani-Schiff (ref. 82) modification to the Baldwin-Lomax (ref. 10) 
turbulence model should improve forebody vortex capture capability when higher aoa solutions 
are computed. Alternative turbulence models should be implemented as they become available 
to improve separation capture. A likely prospect is a model under development patterned off the 
Johnson- King (ref. 11 ) two-dimensional half-equation turbulence model. Along with the improved 
capture of these separated effects, effort should be expended to further develop the topological 
quantification of the interacting fluid mechanisms. Without the development of measures of merit 
for the vortical interactions which take place, the designers’ efforts to control these mechanisms 
will be piecemeal and tentative. 

A final recommendation is to use the TNS codes for other vehicle geometries: This is appro- 
priate both to evaluate the code with alternative data sources, and to begin its realistic application 
as a design tool. Programs already underway include analysis of the transonic flight regimes of 
hypersonic vehicles, and a projected NASA evaluation of the Northrup F-18. This code is being 
seriously considered by several major airframe design companies to assist in the development of 
their Advanced Tactical Fighter (ATF) candidates for the USAF design competition. 
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APPENDIX A 


FLUX JACOBIAN AND VISCOUS COEFFICIENT MATRICES 


The flux jacobians A,B, and C are given by 
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A' 3 u - (7 ~ 1 )AY«’ A'j(7-1) 

K3V — (7 - 1)AY«’ A' 2 (7 — 1) 

A'o + B - A' 3 (7 - 2 )w A' 3 (7 - 1) 

{A'^e/p)-^ 2 ] A'o +7# 
-(7 - l)u>0} 

4 > 2 = 0.5(7 _ l)(^ 2 + V 2 + w 2 ) 

B = K\u + K 2 v + A' 3 u’ 


and, for example to obtain A, 

Ao = it-, Aj (i) A o 
also, the viscous coefficient matrix is given by 


A' 3 - 6 



- 0 

0 

0 

0 

0 


m 2 ] 


q 2^c(/ ? ~ 1 ) 


0 

II 

c 

mzi 


«4 *<(/>“ *) 


0 


m A\ 



«6^c(P _1 ) 

0 


-™51 

W 52 

W53 

”?54 



(A.l) 


(A.2) 


51 


with 


m 2 1 = a.\f)({—u/p) + a 2 S <: {-v/p) + a^b^-w/p) 
m 3i = a 2 b<;(-u/p) + ot 4 b<;(-v/p) + asb^-w/p) 
m 4 i = a 3 b^(-u/p) + a 5 b<;{-v/p) + a^(-w/p) 
m 51 = a\b^{~u 2 / p) + a 2 b((- 2 uv/p) + a 3 b^(- 2 uw / p) 
+ £*4 b^{-v 2 /p) + a 6 b ( (-w 2 / p) + a h b^{- 2 vw / p) 

+ aob^-c/p 2 ) + a 0 ^[(u 2 + v 2 + w 2 )/p) 
m 52 = -m 2 1 - a 0 b^(u/p), m 53 = -m 3] - a 0 b ( {v/p) 

m 54 = -w 4 i - a 0 b^w/p) 

a 0 = jKPr-'id + (y + C s 2 ),«i = Ml(4/3K 2 + C y + Cj] 
a 2 = (A*/3)C*Cy» q 3 = (^/3)C*Cs 
a 4 = p[(i + (4/3 )Cy + Cj], «5 = (A */3)C»C. 
a 6 = a 4C* + Cy + (4/3)Cf] 


(4.3) 


52 


APPENDIX B 


SIMILARITY TRANSFORMATION MATRICES 


The similarity transformation matrix for the Jacobian matrices A, 5, and C is 


T k = 


k x u 


k x v + k z p 
k x w — k y p 

\h<t>V h - l) 


fCy U £ P 

ICyV 

kyW + k x p 

MVb-i) 


k z tl + kyP 
k z v — k x p 

h z w 

[M7(i - 1) 


L +p(k z v - k y w)} +p(k x w - k z u)} +p(k y u - k x v)} 


a 


a 


a(u + k x c) a(u — k x c) 

a(i> + k y c) a(n — k y c) 

a(w + k z c) a(w — k.c) 

a[(4> 2 + c 2 )/( 7 - 1 ) +c6] a{{4> 2 + c 2 )/(~, - 1 ) - c6] J 


(5.1) 


where, k = £, 77 , and ( for A, 5, and C respectively. Also, the inverse transformation matrix is 
given by 

M7 - 


7 1-1 — 
J k ~ 


[k x (l_- 4> 2 /c 2 ) 
—p~ 1 (k z v — fcyte)] 
[k y (l_- <j> 2 /c 2 ) 
—p~ 1 (k x w — A: ; w)] 
\k z {\-4> 2 lc 2 ) 
-p-'ikyU ~ fc,»)] 
iS(4> 7 ~ c6) 
l3{4> 2 + c9) 


. [-kp' 1 

. \-kyP~ 1 
+k z ( 7 - l)nc -2 ] 

/3[k x c - (7 - 1 )u] 

-/3{k x c + (7 + 1 )u] —fS[k y c + (7 — 1 )r] 


[k z p-' 

+M 7 - l)t’C 2 ] 
M 7 - l)t’c -2 

[kxp' 1 

+k z ( 7 - l)nc -2 ] 
t 3 [kyC-( 7 - l)v] 


[- kyP 1 + k x (j-l )wc 2 ] — I‘i(7 — l)c 2 

\kxp ~ 1 + k y ( 7 - l)mc -2 ] -k y ( 7 - l)c -2 

kz(~f — 1 )wc~ 2 —k.(i - l)c~ 2 

i3[k z c - (7 - l)u>] /?( 7-1) 

-/J[Lc + (7 - l)u>] ^(7-1) 

where, 0 = k r u + and, for example, fc, - k x /(k 2 + k 2 4 - k 2 ) 1//2 , etc. 


(5.2) 
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COMPUTATIONAL DOMAIN 


Fig. 2. General-curvilinear coordinate transformation from physical to compu- 
tational space [44]. 
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Fig. 3. Comparison of body-fitted to interfering, non-aligned grids. 
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Fig. 4. Zonal-interface alternatives. 
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Fig. 7. Three-view diagram of the General Dynamics F-16A aircraft [39]. 
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Fig. 8. Examples of: a) unreflexed (or Gothic) and b) reflexed strake planforms. 
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Fig. 9. Surface geometry compared to surface computational grid [40]. 
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Fig. 10. Elliptically generated base-flowfield grid [42]. 
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Fig. 11. Parabolically generated base-flowfield grid [40]. 
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Fig. 12. Front-quarter and rear-quarter views of the surface grid [42]. 
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SYMMETRY PLANE 



Fig. 13. Two-dimensional slice of flowfield grid cross-section normal to stream- 
wise direction [42]. 
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Fig. 14. Planform view of coarse-flowfield grid [42]. 
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Fig. 15. Two-dimensional slice of flowfield cylindrical surface showing airfoil at 
mid-span [42]. 
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Fig. 16. Forebody fuselage cylinder a) on the surface, and b) one station normal 
to the surface [42]. 
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Fig. 17. Chordwise slice of four-zone wing grid at mid-span [44], 
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Fig. 18. Schematic of zonal boundaries on the aircraft surface and in the plane 
containing the wing chordlines (Cuts a, b and c are depicted in Fig. 19). 
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(b) NKC = 7 


Fig. 20. Schematics of zoned circumferential grid slices corresponding to 2 and 
7 base-grid-stations normal to the aircraft surface. 
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Fig. 21. Cell frameworks a) without and b) with clustered interior points. 
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Fig. 23. Clustered grid slice containing airfoil (matches schematic of Fig. 20b). 
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(a) TWO-ZONE GRID SHOWING OVERLAP AT ABCD AND 
EFGH PLANES IN PHYSICAL SPACE 
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(b) GRID POINT DETAIL IN THE OVERLAP REGION IN 
TRANSFORMED SPACE 


Fig. 24. Interface procedure between fine and coarse zones [44]. 
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L2 NORM OF RES 



Fig. 25. Sample convergence rates of the diagonal and block ADI methods 
four-zone solution of a NACA 0012 wing at Moo = 0.826, a = 2° [50]. 



UPPER ZONE 



Fig. 26. Special interface region between zones which straddle the plane 
taining the airfoil chordline. 
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Fig. 27. Cross section of fuselage forebody showing zonal interfaces at plane 
split by strake edge. 


87 



FIRST ZONE 
18 STATION 
FROM INTER- 
FACE PLANE 

INTERFACE 

PLANE 

FIRST ZONE 
17 STATION 
FROM INTER 
FACE PLANE 



Fig. 28. Detail of Fig. 27 showing strake-edge interface between ZONES 17 and 
18. 
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Fig. 29. Classification and depiction of critical points [73]. 
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Fig. 30. A combination of two nodal points of attachment and a saddle point 
[70]. 
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Fig. 32. Behavior of F(y) at large incidence [82]. 
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Fig. 33. Illustration of “streamlines” and “sectional streamlines” [44]. 
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Fig. 34. Coefficient of pressure (C p ) comparisons; Moo = .9, a = 6.15°, 
Re = 4.5 xlO 6 (turbulent solution). 
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Fig. 35. Coefficient of pressure (C p ) comparisons; Mqq = .9, a = 4.12°, 
Re = 4.5 xlO 6 (turbulent solution). 
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Fig. 36 Computed pressure contours on the wing-fuselage planform; M m = 0.9, 
a = 6.15°, Re = 4.5 xlO 6 (turbulent case). 



Fig. 37 Computed Mach contours on the wing-fuselage planform: Moo = 0.9, 
a = 6.15°, Re 4.5 = xlO 6 (turbulent case). 
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Fig. 38. F-16A wind-tunnel oil-flow visualization [39]; Moo = 0.9, a = 6.15°, 
Re = 4.5 xlO°. 
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Fig. 40. Computed surface-particle traces; Mqq = 0.9, a — 6.15°, Re = 4.5 xlO 6 . 


99 



ORIGINAL’ PAGE IS 
OF POOR QUALITY 


Fig. 41. Detail of computed surface-particle traces in strake-wing juncture; 
Moo = 0.9, a = 6.15°, Re = 4.5 xlO 6 . 
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Fig. 42 Computed cross-flow vorticity contours above the strake; Moo = 0.9, 
a = 6.15°, Re = 4.5 xlO 6 (turbulent case). 



Fig. 43 Computed pressure contours above the strake; Moo = 0.9, a = 6.15°, 
Re = 4.5 xlO 6 (turbulent case). 
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Fig. 44 Computed cross-flow velocity vectors above the strake; = 0.9, a = 6.15°, 
Re = 4.5 xlO 6 (turbulent case). 



Fig. 45 Computed cross-flow velocity vectors above the strake; = 0.6, a = 10.0°, 
Re = 8.0 XlO 6 (quasi-laminar case). 
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Fig. 46. Cross-flow topology of vortical interactions above strake. 
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Fig. 47 Computed cross-flow vorticity contours above the strake; Moo = 0.6, 
a = 10.0°, Re = 8.0 xlO 6 (quasi-laminar case). 


Fig. 48 Computed pressure contours above the strake; Moo = 0.6, a = 10.0°, 
Re = 8.0 XlO 6 (quasi-laminar case). 
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Fig. 49 Unrestricted vortex-interaction particle traces, for Moo = 0.6, a = 10.0°, 
Re = 8.0 XlO 6 (quasi-laminar case). 



Fig. 50 F-16A-like water-tunnel flow visualization; a = 15.0°, Re = 1.5 xlO 4 
(Ref. 27). 
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Fig. 51 . Vortex interaction schematic in cross-flow plane at strake-wing juncture. 
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Fig. 52 F-16A water-tunnel flow visualization, for a = 10.0°, Re = 1.5 xlO 4 , 
miniet = 0 (Ref. 21). 



Fig. 53 F-16A water-tunnel flow visualization, for oc = 10.0°, Re = 1.5 xlO 4 , 
riiiniet = 1-0 (Ref. 21). 
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Fig. 54 Unrestricted vortex-interaction particle traces, for Moo = 0*6, oc — 10.0°, 
Re = 1.5 xlO 4 (laminar). 
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Fig. 56. Double-delta configuration water-tunnel flow visualization [91]; 
15.0°, Re = 1.5xl0 4 . 



Fig. 57. F-16A water-tunnel flow visualization [89]; a = 8.0°, Re = 3.0 xlO 4 
per foot. 
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Fig. 58. Locally-restricted surface-particle traces; Moo = 0.6, a. = 10.0°, 
Re = 8.0 XlO 8 (quasi-laminar case). 
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Fig. 59. Critical point analysis of surface-particle traces in Fig. 58. 
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Fig. 60. Detail of Fig. 58 for strake- wing juncture region; Moo = 0.6, a = 10.0°, 
Re = 8.0 xlO 6 (quasi- laminar case) . 
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Fig. 61. Critical point analysis of surface-particle traces on strake in Fig. 60. 
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Fig. 62. Pre-YF-16A wind-tunnel oil-flow visualization; Moo = 0.9, a = 10.71°, 
Re = 5.0 xlO 8 . 
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Fig. 63. Vortical schematic of primary flow properties above the wing; 
Moo = 0.6, a = 10.0°, Re = 8.0 xlO 8 (quasi-laminar case). 



Fig. 64. Locally-restricted surface-particle traces on F16-A geometry; 
Moo = 0.6, ot = 10.0°, Re = 1.5 XlO 4 (laminar case). 
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Fig. 65. Locally-restricted surface-particle traces on strake-wing detail; 
Mqo = 0.6, a — 10.0°, Re = 1.5 xlO 4 (laminar case). 



Fig. 66. Locally-restricted surface-particle traces on F16-A geometry; 
Mqo = 0.3, a = 10.0°, Re = 1.5 xlO 4 (laminar case). 
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